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Summary 

The  project  was  concentrated  on  development  of  new  methodologies  for  decision  mak¬ 
ing  in  uncertain  environment  and  relevant  applications. 

The  first  part  of  the  project  was  focused  on  analytical  and  discrete  optimization  ap¬ 
proaches  for  routing  an  aircraft  in  threat  environment.  The  model  considered  aircraft 
trajectory  in  three-dimensional  space.  Several  threats  were  studied,  including  risk  of 
aircraft  detection  by  radars,  sensors,  and  the  risk  of  being  killed  by  surface  to  air  mis¬ 
siles.  The  problem  of  finding  aircraft  optimal  risk  trajectory  subject  to  a  constraint  on 
the  trajectory  length  was  solved  by  analytical  and  discrete  optimization  approaches. 

The  second  part  of  the  project  resulted  in  general  approach  to  risk  management  for 
the  case  with  uncertainties  in  distributions.  The  risk  of  loss,  damage,  or  failure  was 
measured  by  the  Conditional  Value-at-Risk  (CVaR)  measure.  As  a  function  of  decision 
variables,  CVaR  is  convex,  and  therefore  can  be  efficiently  controlled/optimized  using 
convex  or  linear  programming.  The  methodology  was  tested  on  two  Weapon-Target 
Assignment  (WTA)  problems.  The  total  cost  of  a  mission  was  minimized,  while  sat¬ 
isfying  the  operational  constraints  and  ensuring  destruction  of  targets  with  high  prob¬ 
ability.  The  risk  of  a  failure  of  the  mission  is  controlled  by  CVaR  constraints.  The 
case  studies  showed  that  there  were  significant  qualitative  and  quantitative  differences 
in  solutions  of  deterministic  and  stochastic  WTA  problems. 

The  third  part  of  the  project  studied  the  Multiple  Traveling  Salesmen  Problem 
(Multiple-TSP)  in  several  variations.  The  research  was  focused  on  MIN-MAX  2-TSP, 
which  cannot  be  solved  by  standard  methods.  The  relation  between  this  class  of  prob¬ 
lems  and  a  subclass  of  the  self-dual  monotonic  Boolean  functions  was  established.  This 
resulted  in  new  efficient  optimization  algorithms. 
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1.  Introduction 

The  class  of  military  and  civil  engineering  applications  dealing  with  opti¬ 
mal  trajectory  generation  for  space,  air,  naval  and  land  vehicles  is  very  broad. 
It  addresses  several  types  of  problems  with  various  objectives,  constraints  on 
resources  and  control  limitations,  for  instance, 

■  Minimizing  risk  of  aircraft  detection  by  radars,  sensors  or  surface  air 
missiles  (SAM)  [5, 19, 22] 

■  Minimizing  risk  of  submarine  detection  by  sensors  [21] 

■  Minimizing  cumulative  radiation  damage  in  passing  through  a  contami¬ 
nated  area 

m  Finding  optimal  trajectories  for  multiple  aircraft  avoiding  collisions  [15] 

■  Maximizing  the  probability  of  target  detecting  by  a  searcher  [1,  3,  9, 12, 
13,16,17,20] 

■  Minimizing  propellant  consumption  by  a  spacecraft  in  interplanetary 
and  orbit  transfers  [4] 

■  Minimizing  a  weighted  sum  of  fuel  cost  and  time  cost  for  a  commercial 
plane 

■  Minimizing  energy  for  a  mobile  robot  on  terrains 

We  are  interested  in  developing  efficient  optimization  approaches  capable 
of  solving  a  broad  class  of  applications  related  to  trajectory  optimization.  This 
chapter,  being  the  first  step  in  accomplishing  this  task,  is  primarily  focused 
on  optimal  path  planning  for  an  aircraft  in  a  threat  environment.  The  threat 
is  associated  with  the  risk  of  aircraft  detection  by  radars,  sensors  or  SAMs. 
The  chapter  develops  analytical  and  discrete  optimization  approaches  to  op¬ 
timal  trajectory  generation  that  minimize  the  risk  of  aircraft  detection  with: 
1)  variable  aircraft  Radar  Cross-Section  (RCS);  2)  different  types  of  detecting 
installations;  3)  arbitrary  number  of  detecting  installations;  4)  constraint  on 
trajectory  length;  and  suggests  efficient  algorithms  for  solving  the  formulated 
risk  minimization  problem. 

Optimal  trajectory  generation  is  a  fundamental  requirement  for  military  air¬ 
craft  flight  management  systems.  These  systems  are  required  to  take  advantage 
of  all  available  information  in  order  to  perform  integrated  task  processing,  re¬ 
duce  pilot  workloads  and  provide  updates  at  regular  time  intervals  sufficient 
for  threat  avoidance  [19].  A  model  for  optimal  routing  an  aircraft  in  a  threat 
environment  is  developed  based  on  specified  mission  objectives,  available  re¬ 
sources  (fuel  capacity),  limitations  on  aircraft  control  while  minimizing  risk 


6 


exposure.  In  general,  it  addresses  uncertainty  and  dynamics  inherent  to  opti¬ 
mal  path  planning  and  makes  idealizing  assumptions  with  respect  to  geomet¬ 
rical  and  physical  properties  of  an  aircraft  and  threat  environment.  Despite 
numerous  studies  in  this  area,  only  a  few  considered  risk  optimization  prob¬ 
lems  with  technological  constraints.  Zabarankin  et  al.  [22]  suggested  analyt¬ 
ical  and  discrete  optimization  approaches  for  optimal  risk  path  generation  in 
two-dimensional  (2D)  space  with  constant  RCS,  arbitrary  number  of  sensors 
and  a  constraint  on  path  length. 

This  chapter  develops  a  3D  model  for  minimizing  risk  of  aircraft  detection 
by  radars,  sensors  or  SAMs  with  variable  RCS.  A  sensor  is  considered  to  be  an 
antenna  capable  of  receiving  an  isotropically  radiated  signal  from  the  aircraft, 
while  a  radar  is  assumed  to  be  an  antenna  capable  of  transmitting  a  signal  and 
receiving  the  signal  reflected  off  of  the  aircraft.  The  model  is  deterministic  and 
static,  since  it  assumes  no  uncertainty  in  aircraft  detection  and  radar  locations 
and  considers  neither  aircraft  kinematics  equations  nor  parameters  for  aircraft 
cohtrol  during  a  flight.  The  risk  of  detection  is  assumed  to  be  independent  on 
aircraft  speed.  This  model  extends  the  2D  risk  minimization  problem  [22]  of 
aircraft  detection  by  sensors  to 

■  3D  space 

■  Variable  RCS  —  an  aircraft  is  considered  to  bean  axisymmetrical  ellip¬ 
soid  with  the  axis  of  ellipsoid  symmetry  determining  direction  of  aircraft 
trajectory 

■  Risk  of  detection  to  be  proportional  to  the  aircraft’s  RCS  and  reciprocal 
to  the  rfi^-power  of  the  distance  between  the  aircraft  and  a  particular 
detecting  installation,  where  n  =  2  corresponds  to  a  passive  listener  or 
sensor,  and  n  =  4  corresponds  to  an  active  listener  or  radar 

The  purpose  of  this  .simplified  model  is  analyzing  the  impact  of  variable 
RCS  on  the  3D  geometry  of  optimal  trajectories  subject  to  a  constraint  on  tra¬ 
jectory  length  and  evaluating  performance  of  the  developed  discrete  optimiza¬ 
tion  approach  with  respect  to  running  time  and  accuracy.  Verified  optimization 
techniques  will  be  applied  in  optimal  path  planning  with  actual-tabulated  RCS. 

We  developed  analytical  and  discrete  optimization  approaches  for  solving 
formulated  trajectory  optimization  problem  with  a  constraint  on  trajectory  length 
and  arbitrary  number  of  detecting  installations  (sensors  or  radars).  Through 
techniques  of  Calculus  of  Variations,  the  necessary  optimality  conditions  for  a 
solution  to  the  risk  minimization  problem  were  reduced  to  a  nonlinear  vectorial 
differential  equation.  In  the  case  of  a  single  radar,  we  obtained  an  analytical 
solution  to  this  equation  expressed  by  a  quadrature.  Analytical  solutions  are 
intended  for  conceptual  understanding  and  analyzing  the  impact  of  variation 
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in  RCS  on  the  geometry  of  optimal  trajectories  and  testing  performance  of  the 
developed  discrete  optimization  approach  in  the  case  of  a  single  radar  in  2D 
and  3D  spaces.  Although  we  have  made  significant  progress  in  the  develop¬ 
ment  of  &e  analytical  approach,  finding  an  analytical  solution  to  the  vectorial 
differential  equation  in  the  case  of  an  arbitrary  number  of  installations  is  still 
an  open  issue.  This  is  one  of  the  main  reasons  for  addressing  development  of 
discrete  optimization  approaches. 

Several  discrete  optimization  approaches  are  available  for  numerical  solving 
proposed  risk  minimization  model.  All  these  approaches  may  tentatively  be 
divided  into  three  major  categories: 

■  Gradient-based  algorithms 

m  Dynamic  programming 

■  Network  Flow  (NF)  optimization 

Efficiency  of  discrete  optimization  approaches  in  optimal  risk  path  planning 
essentially  depends  on  type  of  risk  functionals,  technological  constraints,  and 
a  scheme  of  trajectory  approximation  (see,  for  instance,  [19]  for  discussions 
of  these  issues).  Gradient-based  algorithms  are  very  efficient  when  the  risk  of 
detection  is  determined  by  smooth  analytical  functionals.  However,  while  dy¬ 
namic  programming  and  NF  optimization  are  global  optimization  approaches, 
gradient-based  algorithms  most  likely  find  only  locally  optimal  solutions  in 
die  case  when  risk  functionals  are  nonconvex.  Many  of  the  previous  studies  on 
trajectory  generation  for  military  aircraft  are  concentrated  on  feasible  direction 
algorithms  and  dynamic  programming  [5].  These  methods  tend  to  be  compu¬ 
tationally  intense  and,  therefore,  are  not  well  suited  for  onboard  applications. 
To  improve  computation  time,  John  and  Moore  [19]  used  simple  analytical  risk 
functions.  Based  on  such  an  approach,  they  developed  lateral  and  vertical  al¬ 
gorithms  to  optimize  flight  trajectory  with  respect  to  time,  fuel,  aircraft  final 
position,  and  risk  exposure.  Nevertheless,  these  algorithms  are  not  intended  for 
solving  optimization  problems  with  technological  constraints,  such  as  a  con¬ 
straint  on  the  trajectory  length.  Zabarankin  et  al.  [22]  demonstrated  efficiency 
of  NF  optimization  approach  in  solving  risk  minimization  problem  with  a  con¬ 
straint  on  trajectory  length  and  arbitrary  number  of  sensors  in  2D  space.  The 
main  advantages  of  using  NF  optimization  approach  are 

■  Among  all  feasible  approximated  trajectories  in  a  considered  network, 
NF  approach  finds  a  globally  optimal  one. 

■  Complexity  of  NF  algorithms  is  independent  on  number  of  detecting 
installations. 

■  It  can  easily  be  applied  for  the  case  with  actual  nonsmooth  RCS. 
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Despite  these  advantages,  the  complexity  of  NF  algorithms  substantially 
depends  on  the  coarseness  of  a  network,  in  particular,  on  the  number  of  arcs. 
Consequently,  precision  for  an  optimal  solution  should  reasonably  be  specified. 
Recently,  Tsitsiklis  [18]  and  Polymenakos  et  al.  [14]  suggested  Dijkstra-like 
and  corfecting-like  methods  for  efficient  solving  a  continuous-space  shortest 
path  problem  in  2D  plane.  In  this  case,  finding  a  globally  optimal  trajectory 
employs  discretization  of  Hamilton- Jacobi  equation  [18],  which  turns  out  to  be 
an  efficient  synthesis  of  analytical  and  discrete  optimization  techniques.  This 
supports  the  philosophy  that  using  analytical  properties  of  objective  functions 
in  NF  optimization  leads  to  more  efficient  algorithms.  Since  our  goal  is  gen¬ 
erating  globally  optimal  trajectories  and  on  the  next  step  applying  developed 
optimization  approach  in  optimal  path  planning  with  actual-nonsmooth  RCS 
(in  this  case  utilizing  analytical  properties  of  risk  functionals  is  limited),  we 
considered  NF  optimization  approach. 

We  approximated  an  admissible  domain  for  aircraft  trajectory  by  a  3D  net¬ 
work  with  a  flexible  structure  and  presented  aircraft  trajectory  by  a  path  in 
this  network.  NF  optimization  approach  reduced  optimal  risk  path  generation 
with  a  constraint  on  trajectory  length  to  the  Constrained  Shortest  Path  Prob¬ 
lem  (CSPP).  Development  of  efficient  network  structures  with  relatively  small 
numbers  of  arcs  and  nodes  while  preserving  flexibility  for  trajectory  approxi¬ 
mation  is  one  of  the  key  issues  in  reduction  of  approach  computational  time. 

To  solve  the  CSPP  in  2D  and  3D  cases,  we  used  the  Label  Setting  Algorithm 
(LSA)  with  a  preprocessing  procedure  [7,  8]  and  network  structure  smooth¬ 
ing.  The  efficiency  of  the  discrete  optimization  approach  is  demonstrated  by 
several  numerical  examples  with  various  ellipsoid  shapes,  constraints  on  tra¬ 
jectory  length  in  the  cases  of  one,  two  and  three  radars.  For  the  case  with  a 
single  radar,  we  compared  analytical  and  numerical  solutions  and  found  that 
solutions  coincide  with  high  precision  in  2D  case  and  are  very  close  in  3D 
case.  The  fact  that  discrete  trajectories  are  closer  to  corresponding  analytical 
ones  in  2D  case  can  be  explained  by  different  flexibility  of  2D  and  3D  network 
structures  in  trajectory  approximation.  LSA  running  time  in  all  2D  testing  ex¬ 
amples  is  only  several  seconds,  indicating  that  this  NF  algorithm  is  fast  enough 
for  use  in  online  applications  with  a  relatively  small  number  of  arcs  in  a  graph. 
However,  it  is  also  known  that  the  CSPP  is  an  NP-hard  problem  and,  conse¬ 
quently,  no  exact  polynomial  algorithms  should  be  expected.  Numerical  tests 
in  a  3D  case  reveal  that  LSA  running  time  strongly  depends  on  the  shape  of 
ellipsoid.  This  phenomenon  has  been  analyzed  from  optimization  perspective 
and  an  improvement  for  preprocessing  procedure  has  been  suggested. 

The  chapter  is  organized  as  follows:  section  2  develops  a  3D  model  for 
trajectory  optimization  with  variable  RCS  subject  to  a  constraint  on  trajec¬ 
tory  length;  section  3  derives  the  vectorial  differential  equation  for  finding  the 
optimal  trajectory  in  a  general  case  and  obtains  analytical  solution  to  this  equa- 
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tion  in  the  case  of  a  single  radar;  section  4  reduces  optimal  path  planning  to 
the  CSPP  and  presents  the  LSA  with  preprocessing  procedure  and  smooth¬ 
ing  condition;  section  5  conducts  numerical  experiments  with  various  ellipsoid 
shapes  and  constraints  on  trajectory  length  in  the  cases  of  one,  two  and  three 
radars;  section  6  analyzes  results  of  numerical  experiments  from  optimization 
and  variable  RCS  perspectives;  section  7  discusses  main  analytical  and  nu¬ 
merical  results  and  concludes  the  chapter;  the  appendix  considers  necessary 
optimal  conditions  for  calculus  of  variations  problem  with  a  nonholonomic 
constraint  and  movable  end  point. 


2.  Model  Development 

This  section  develops  a  three-dimensional  (3D)  model  for  minimizing  the 
risk  of  aircraft  detection  by  a  network  of  active  or  passive  installations  (radars, 
sensors)  with  variable  aircraft  RCS. 

Suppose  an  aircraft  must  fly  from  point  A,  {xaiVAiZa),  to  point  B, 
{xbiVbi^b)^  in  3D  space  trying  to  minimize  the  cost  of  detection  from  N 
radars  located  in  the  area  of  interest.  We  model  the  aircraft  by  an  axisymmetri- 
cal  ellipsoid  with  the  axes’  lengths  a,  h  and  6.  The  axis  with  length  o  is  the  axis 
of  ellipsoid  symmetry,  which  orients  a  direction  of  aircraft  trajectory.  Ellipsoid 
shape  is  defined  by  parameter  k  =  b/ a.  Cases  of  k  =  1,  k  <  1  and  /c  >  1 
correspond  to  sphere,  elongated  and  compressed  ellipsoids,  respectively,  see 
Figure  1.1. 


Figure  1.1.  Ellipsoid  shape  is  defined  by  parameter  k  =  h/a. 

Let  vectors  r  =  {x,y,z)  and  qi  =  {ai,bi,Ci),  i  =  1,N,  determine  po¬ 
sition  of  ellipsoid  geometrical  center  and  position  of  the  t***  radar,  respec¬ 
tively.  A  trajectory  of  the  ellipsoid’s  center  is  assumed  to  be  a  path  of  the 
aircraft.  We  define  a  trajectory  as  a  function  of  its  current  length  s,  i.e.  r  = 


10 


r(s)  =  (x(s),y{s),z{s)).  Such  a  parameterization  is  also  known  to  be  the 
natural  definition  of  a  curve.  Vector  f(s)  =  ^r(s)  =  (a:(s),y(s),i(s))  de¬ 
termines  a  direction  of  aircraft  trajectory  that  coincides  with  the  axis  of  el¬ 
lipsoid  symmetry.  Since  (ds)^  =  {dx)^  +  {dyf'  H-  {dz)"^ ,  vector  f (s)  must 
satisfy  condition  =  1.  The  length  of  vector  ri(s)  = 

r(s)  -  q*  =  {x-ai,y-  bi,z-  Ci),  denoted  by  l|ri(s)|l,  defines  the  dis¬ 
tance  from  the  aircraft  to  the  installation  (see  Figure  1.1),  i.e.  ||ri(s)||  = 


X 


Figure  L2.  3D  model  for  optimal  path  planning  in  a  threat  environment. 

The  RCS  of  the  aircraft  exposed  to  the  i***  radar  at  point  {x,  y,  z)  is  propor¬ 
tional  to  the  area  of  the  ellipsoid’s  projection  to  the  plane  orthogonal  to  vector 
r< 

RCSj  =  aiSi, 

where  the  constant  coefficient  ai  depends  on  the  radar’s  technical  character¬ 
istics  such  as  the  maximum  detection  range,  the  minimum  detectable  signal, 
the  transmitting  power  of  the  antenna,  the  antenna  gain  and  the  wavelength  of 
radar  energy. 

The  magnitude  of  ellipsoid  projection  area  is  given  by  the  formula  Si  = 
TT  b  y/a^  sin^  -f-  6^  cos2  Ou  where  Oi  is  the  angle  between  vectors  and  r. 
Based  on  relation  cos  9i  =  and  using  notation  k  =  b/a,  the  formula  for 
RCSt  is  identically  rewritten  as 
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The  purpose  of  presenting  RCSi  in  the  form  of  (1)  is  the  following.  Since 
the  aircraft  has  a  limited  size,  we  assume  the  value  ^/aF+~^  (’’diameter”  of 
cross-section)  to  be  constant  for  all  a  and  b  and,  hence,  the  form  of  the  ellipsoid 
to  be  defined  by  ratio  6/a,  i.e.  parameter  k,  only.  For  instance,  the  case  of 
K  =  6/a  =  0  corresponds  to  an  infinitely  thin  needle  with  the  length  of  a, 
while  the  case  of/^  =  6/a  — ^  oo  corresponds  to  an  infinitely  thin  disk  with 
the  radius  of  6.  Note  that  the  cross-section  of  the  infinitely  thin  needle  always 
equals  zero,  whereas  the  cross-section  of  the  infinitely  thin  disk  is  reduced  to 
aiTta^l  cos  6i\,  which  is  zero  only  when  =  f . 

The  risk  function  (also  referred  to  as  cost  function)  for  detection  of  the  air¬ 
craft  by  the  radar  is  proportional  to  the  RCS  and  reciprocal  to  the 
,  power  of  the  distance  between  the  aircraft  and  the  radar,  ||ri||’^  (cases  of 
n  =  2  and  n  =  4  correspond  to  sensor  and  radar,  respectively),  namely, 


RCSi  fa‘^  +  b‘^\ 


2k 


1  +  (/c2  -  1) 


+  K"* 


Since  the  value  of  y/cfi  +  is  assumed  to  be  constant,  product  tt  ^  can 
be  omitted  for  simplicity  and  the  risk  function  for  detection  of  the  aircraft  by 
the  i***  installation  is  reexpressed  with  normalized  coefficient  Wi 


C(ri,f) 


2k Wi  /ikilP  +  1)  (ti  •  f ) 

1+^  ||ri|h+i 


(2) 


/  N  N 

where  Wi  =Oi  j  and  ^  lOj  =  1. 

/  t=l  i=l  .  /  N  t. 

We  assume  the  risk  of  detection  from  iV  radars  at  point  r  =  (x,  y,  z)  to  be 
the  sum  of  risk  functions  (2)  for  all  i  =  1,  iV 


N 


L(r,r)  =  ^C(ri,r)  = 

t=l  1=1 


2«  i  711^11^  + 


|r,|l”« 


(3) 


The  total  risk  of  detection  is  the  integral  of  (3)  along  aircraft  trajectory  with 
length  I,  i.e. 

l 

^(r,r)  =  J  L{r{s),T{s))  ds. 

0 


The  risk  minimization  problem  is  finding  a  trajectory 

p  =  r{s)  =  {x{s),y{s),z{s))  ,0<s<l, 


(4) 
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from  point  A  to  point  B,  having  coordinates  =  (xa,  Pa,  ^a)  and  = 
(xb,  VBi  zb)j  respectively,  with  the  minimal  risk  of  detection  subject  to  a  con¬ 
straint  on  trajectory  length 

min  ^(r.f) 

(5, 

r(0)  =  ta,  r(0  =  fb, 
l<l*. 

The  form  of  risk  functional  (4)  implies  that  either  the  risk  is  independent 
on  aircraft  speed,  or  aircraft  speed  is  always  a  unit.  Under  assumption  of  unit 

speed,  $  becomes  time  variable  t,  total  length  I  becomes  total  time  T,  and  (5) 

T 

is  viewed  as  a  problem  of  optimal  control  with  .F(r,  v)  =  /  L  (r(f),  v(f))  dt, 

0 

with  V  =  r.  Whatever  the  interpretation  of  (5)  is,  analysis  techniques  are  the 
same.  To  solve  problem  (5),  calculus  of  variations  and  network  flow  optimiza¬ 
tion  approaches  are  addressed. 

We  want  to  mention  the  2D  dynamic  model  for  minimizing  the  risk  of  sub¬ 
marine  detection  by  a  network  of  sensors  [21].  In  that  model,  the  risk  func¬ 
tional  considers  different  directions  of  the  power  radiation  and  variable  speed 
of  a  submarine 


f  q  0-~  K2COs(2di(t)))  , 
0 


where  ri(i)  and  &i(t)  mean  exactly  the  same  as  in  aircraft  detection  model 
(see  Figure  1.2),  i.e.  ri(t)  is  the  distance  between  the  submarine  and  the 
i***  sensor  and  0i(t)  is  the  corresponding  angle,  which  now  both  depend  on 
time  t;  rj  denotes  submarine  speed  at  time  moment  t;  Si  is  the  sensitivity 
coefficient  associated  with  the  i***  sensor;  ki  and  K2  are  parameters,  which 
correspond  to  doubling  the  radiated  power  relative  to  small  speeds  (usually, 
Ki  =  0.0003  <<  1)  and  adjustment  of  power  radiation  in  different  directions 
(k2  <  1;  when  K2  >  0,  power  tends  to  be  radiated  most  strongly  broadside; 
K2  is  arbitrarily  set  to  0.5  [21]),  respectively.  The  term  1  —  /C2  cos(20i(f))  may 
be  considered  as  submarine  cross-section  in  2D  space.  The  Optimal  Control 
approach,  suggested  to  solve  the  model,  starts  with  some  feasible  trajectory 
provided  by  an  observer  and  transforms  it  to  locally  optimal  one  by  steepest- 
descent  technique. 


3.  Calculus  of  Variations  Approach 

This  section  presents  a  vectorial  differential  equation  for  solving  the  risk 
minimization  problem  (5)  and  obtains  an  analytical  solution  to  this  problem 
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in  the  case  of  a  single  radar  or  sensor.  The  vectorial  differential  equation  is 
the  reduction  of  necessary  conditions  for  an  extremal  minimizing  a  functional 
with  a  nonholonomic  constraint  and  movable  end  point.  Through  techniques 
of  Calculus  of  Variations,  the  appendix  derives  this  equation  in  a  general  case. 

Introducing  a  new  function 

9i{ri,  r,  Xl)  = - ^ ^  (6) 

and  using  notation  gi  =  we  formulate  necessary  conditions  for  an  opti¬ 
mal  trajectory  satisfying  (5). 

Theorem  3.1.  (vectorial  differential  equation).  An  optimal  solution  to  opti¬ 
mization  problem  (5)  should  necessarily  satisfy  the  following  vectorial  differ¬ 
ential  equation 


with  boundary  conditions 

r(0)  =  fa,  r(Z)  =  r^,  l<  I*, 
and  nonholonomic  constraint 

</){r)  =  —  1  =  0. 


(7) 

(8) 

(9) 


Proof  The  problem  (5)  is  a  particular  case  of  the  problem  (1  .A.  1)-(1  A.5)  con¬ 
sidered  in  the  appendix.  In  the  case  of  (9)  we  have  =  2f  and  consequently 


dd>\  _ 


=  r.  Substituting  the  last  equality  into  the  general  vectorial 
differential  equation  (1.A.10)  derived  in  the  appendix,  we  obtain  the  vectorial 
differential  equation  for  determining  an  optimal  trajectory  r(s),  0<  s  <1, 


dL 

dr 


=  0. 


(10) 


Introducing  a  new  constant  Ax,  by  relation  cx,  =  f^^Ax  and  using  notations 
9i  =  9(ri,  r,  Ax),  gi  =  £g  (r,,  r,  Ax),  where  function  g{ri,  r,  Ax)  is  defined 
by  (6),  we  verify  that  the  following  relations  hold  for  (3) 


L-f 


dL 

dr 


+  cl  = 


2k 
1  -f 


N 

t=l 
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dL  ddL  _  2k  a  .  f  Vj  \ 
dr  ds  dr  1  +  k^  ^  \ri  •  rj  ’ 

which  reduce  (10)  to  equation  (7). 

Note  equation  (7)  and  constraint  (9)  are  dependent  in  the  sense  that  the  scalar 

N 

product  of  (7)  with  f  is  reduced  to  ((l  —  r^)  in  —  which 

i=l 

becomes  identity  if  (9)  is  satisfied.  lj 


Remark:  equation  (7)  may  be  presented  in  different  forms.  Using  relation  (5), 
we  have  ^-r=  and  -r(ri  •  f)  =  f  x  [ri  x  r].  Consequently,  (7) 

becomes 


^  /^r  X  [ri  X  f]  .  ^ 

i=l  ^  1=1 


r,  K  l(f  s.) 


0. 


Equation  (7)  may  also  be  presented  in  a  matrix  form.  Denoting 


Gjvxi  = 


WlQl  \ 

W2g2 

,  e;\rxi  = 

/  1  \ 
1 

WN9N  j 

J 

) 


RzxN 


■fe' 


r2 

r2  •  r’ 


rN  \ 
rjv-f/  ’ 


(7)  is  rewritten  as 

ilG=^(f(e-^.G)). 

Choice  of  a  form  for  (7)  is  just  a  matter  of  convenience  for  conducting  an¬ 
alytical  manipulations  or  numerical  analysis.  Differential  equation  (7)  may  be 
solved  numerically  by  an  appropriate  gradient-based  algorithm,  however,  in 
this  case  we  are  not  guaranteed  to  obtain  a  globally  optimal  solution. 


Deriving  an  analytical  solution  to  (7)-(9)  with  an  arbitrary  number  of  radars 
or  sensors  is  reduced  to  finding  the  second  integral  for  equation  (7)  (the  first 
one  is  nonholonomic  constraint  (9)),  which  still  remains  an  open  issue.  The 
next  theorem  shows  how  the  second  integral  and  a  corresponding  analytical 
solution  are  found  in  the  case  of  a  single  radar  or  sensor. 


Theorem  3.2.  (the  case  of  a  single  radar  or  sensor).  In  the  case  of  a  single 
radar  or  sensor,  located  at  the  origin  of  the  system  of  coordinates,  i.e.  point 
(0,0.0). 
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(1)  the  optimal  trajectory  is  a  planar  curve  in  3D  space,  where  the  trajec¬ 
tory's  plane  is  determined  by  the  origin  of  the  system  of  coordinates  arui  the 
starting  and  finishing  trajectory’s  points,  i.e.  by  (0, 0, 0),  and  tb  (the  equa¬ 
tion  of  the  plane  is  given  by  [r^  x  r^]  •  r  =  0^; 

(2)  introducing  a  polar  system  coordinates  {p,  ip)  in  the  trajectory’s  plane, 
vectorial  differential  equation  (7)  with  (8)  and  (9)  is  reduced  to  a  nonlinear 
first-order  differential  equation  with  respect  to  function  p  =  p{ip) 

- .  ^  - +  .  . -  =  C,  (11) 

/?”-2^/c2(p'^)2  +  p2  y(p^)2  +  p2 

with  boundary  conditions 

p{fpA)  =  PA,  P{'4’b)  =  PB,  (12) 

defining  points  A  and  B  in  the  polar  system  (p,  ip),  and  a  constraint  on  trajec¬ 
tory  length 

‘‘I’B 

f  i/^p)^diP  =  h.  (13) 

^l>A 

Proof.  Since  an  analytical  solution  to  (7)  is  derived  in  the  case  of  a  single 
radar,  without  loss  of  generality,  we  assume  that  the  radar  is  located  at  the 
origin  of  the  system  of  coordinates,  that  is,  (ai,  6i,  ci)  =  (0, 0, 0),  and  ri  = 
r.  Functions  L  (r,  f),  g  (r,  r,  Ax,)  and  equation  (7)  in  this  case  are  presented, 
respectively, 

„  2k  v'lWP  +  (K»-l)(r-f)’= 

Mr,r)-i_j^^2  ||r||n+i 


g  (r,  f ,  Ai)  = - /  . .  + 

l|r||"-Vl|rl|2  +  (/c2-l)(r.f)2 

Producing  vectorial  product  of  the  last  equation  with  vector  r,  we  obtain 

|(l■•x^|s)  =  0, 

which  is  equivalent  to  having  the  first  integral 

[r  xr]p  =  C, 


(14) 
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where  C  =  {C\,C2,  C3)  is  a  constant  vector.  Since  (r  •  [r  x  f])  =  0  and 
g  (r,  t,XL)^  0.  the  scalar  product  of  (14)  with  r  leads  to 

C  •  r  =  0, 

which  is  the  equation  of  a  plane  going  through  the  origin  of  the  system  of 
coordinates.  It  means  that  an  optimal  trajectory  is  a  planar  curve  in  a  3D 
space,  i.e.  all  its  points  form  a  single  plane  in  a  3D  space  (lie  within  the  same 
plane).  Since  boundary  points  A  and  B  must  also  belong  to  trajectory’s  plane, 
i.e.  vectors  ta.  tb  must  satisfy  equation  C  •  r  =  0,  vector  C  is  parallel  to 
[fa  X  fb],  and  the  explicit  expression  for  the  trajectory’s  plane  is  given  by 

[fa  X  fb]  •  r  =  0,  (15) 


or 

hxX  +  hyy  +  hzZ  —  Q., 

where  {hx,  hy,hz)  are  the  components  of  vector  [fa  x  fb] 

{hx,  hy,hz)  =  iVAZB  -  VBZA,  ZAXB  “  ZbXa,  xaVb  “  xbVa)- 


The  next  step  is  parameterizing  3D  plane  (15)  by  a  2D  polar  system  of  coor¬ 
dinates  (p,  Ip).  A  point  with  coordinates  (a;(p,  'ip),  y{p,  ip),z{p,  ip))  should  sat¬ 
isfy  (15)  identically.  Suppose  the  origin  of  the  polar  system  (p,  ip)  coincides 
with  the  origin  of  the  original  3D  system  of  coordinates,  i.e.  point  (0, 0, 0). 
Let  Ip  he  a  counterclockwise  angle  producing  left-handed  screw  with  the  vec¬ 
tor  [fa  X  Fb]  and  counted  from  the  upper  side  of  the  plane  xy.  Introducing 
notations 

^ — gg .  COS  l3  ” 


cos  a  = 


sin  a  = 


hy _ 

v/SI+^’ 


sm 


y/hl+hl+h^,  ’ 

P-  /4+hl+h!i' 


coordinates  (»,  y,  z)  of  points  identically  satisfying  (15)  are  determined  by  the 
following  relations 


x{p,ip)  =  p  (sin  a  cosip  —  cos  a  cos/3  sin  ip) , 
y{p,  Ip)  =  — p  (cos  a  cos  ^  4-  sin  a  cos  /3  sin  ip) , 
z{p,ip)  =  p  sin /3  sin^. 


Based  on  these  relations,  we  have 

[f  X  r]  =  —p^ip 


2.;.  [faxfb] 


||[fa  X  fbIH’ 

and,  consequently,  using  the  last  formula,  (14)  is  reduced  to  the  following 

p^ipg  =  C, 


scalar  equation 


(16) 
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where  C  is  unknown  constant  scalar  value.  Since  =  ||r|p  =  ,  function 

g  (r,  f ,  Xl)  is  rewritten  as 

9  (P,  P,  W  =  „  A  .  2  i'v2 

P”a/1  +  (/t2  -  l)f/ 

With  this  relation,  equation  (16)  and  constraint  (9),  expressed  in  terms  of 
{p,  -ip)  as  p^  +  p^iP'  =  1,  determine  a  system  of  differential  equations  for 
finding  optimal  p{s)  and  'ij){s) 


\  - -  ^  +  Xlp^  ]  tp  =  C, 

with  boundary  conditions  p(0)  =  pa,  ^(0)  =  ‘ipA, 
where  {pA,  V’a)  and  {ps,  iPb)  are  given  by 


p'  +  pV  =  i, 

p{i*)  =  PB, 


PA  =  i’A  =  arccos  . 

PB  =  llrfill,  V-B  =  aiccoa  . 


(17) 


V'B. 


Let  p'^  =  Using  relation  p  =  p'^ip  with  the  second  equation  of  (17)  we 
present  p  and  -0  as 


Substitution  of  the  last  formulas  into  the  first  equation  of  (17)  eliminates 
variable  s  from  the  system  (17)  and  reduces  it  to  the  nonlinear  first-order  dif¬ 
ferential  equation  (11)  determining  p  as  a  function  of  ip  with  boundary  condi¬ 
tions  (12).  Since  variable  s  was  eliminated  from  (17),  the  second  equation  of 
(17)  is  satisfied  identically,  and,  thus,  a  constraint  on  trajectory  length  should 
be  included  in  the  form  of  (13).  .  . 

Note  it  does  not  matter  what  sign,  plus  or  minus,  we  choose  for  ip  in  ip  = 
■£  ^  — ,  since  we  always  can  change  the  sign  of  the  constant  C  in  the 

ip',i,P+p^ 

right-hand  side  of  equation  (14)  and  denote  it  by  a  new  constant.  □ 


Discussion  of  necesssiry  and  sufficient  conditions  for  a  mini¬ 
mum.  Equation  (1 1)  (or  system  (17))  is  only  the  necessary  condition  for  a 
trajectory  to  be  just  an  extremal,  since  (11)  finds  trajectories  minimizing  the 
functional  (4)  under  given  conditions  as  well  as  maximizing  it.  A  sufficient 
condition  for  a  solution  to  system  (11)  to  be  an  optimal  trajectory  (p*,  •0*),  i.e. 
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to  minimize  the  functional,  requires  the  second  variation  of  the  functional  at 
tA*)  to  be  greater  than  or  equal  to  zero  (see  the  appendix). 

In  the  case  of  a  single  radar,  the  risk  functional  with  the  relaxed  constraint 
f?  +  =  1  is  presented  in  polar  coordinates  by 


i 


0 


Lip,  pA) 


2k  +  -1)p^  ^  ^ 

1  +  p" 


(/?2  +  p2^2  _ 


Where  A*  =  ^  +  A^J  =  5*-  Assuming  the  con¬ 

straint  on  the  length  of  a  trajectory  to  be  active,  i.e.  I  =  Z*,  the  second  variation 
of  the  functional  at  (p»,  ■0*)  is  defined 


=  s{L;,(5p)'‘+L;i,{Sf,f+L';i^iH?+2L;fSpp 

+  2LrSp  + 

=  /  ((i«  -  (w* + 

L"  .,6p  i}  +  2L",6p5i;)\  ds 

r  /  {P{5pf  +  Q{Sp)^  +  g*pi{Hf  +  ^9*pA*^P  ds 
/f[(p-4gAi)  {Sp?  +  QiSp? 

+g*  (^ip^Sp  -H  p*#)  ^  ds, 
where  P  and  Q  are  given  by 


^  -w-f  I 

+ 

_  2k 

l+K^ 

=  2k 

1+7? 


p=-+ 


n  d  ( 


Q  = - - t+9- 

p”(l -1- («2  _  l)p2)2 

Since  the  extremal  (p*,  V’*)  satisfies  (17),  we  can  use  (17)  to  rearrange  P ,  Q 
and  the  other  terms  in  toe  integral  of  and,  thus,  obtain  different  equivalent 
expressions  for  S'^P.  However,  verification  of  toe  condition  S^P  >  0  for  all 
<5p,  6p  and  even  in  this  particular  case,  is  not  a  trivial  task. 
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We  confine  ourselves  here  only  to  verification  of  necessary  conditions  for 
a  minimum.  For  an  extremal  to  minimize  the  functional,  it  should 

necessarily  satisfy  the  Legendre  conditions 


lp=p*,7p=lp* 


^pi> 

L"  ■  L- 

f»P  p=p*,'^ip* 


which  are  reduced  to  verification  of  Q  >  0  and  g*  >  0.  In  the  case  of  A:  > 
1,  condition  Aj,  >  0  guarantees  satisfaction  of  both  Q  >  0  and  g*  >  0, 
which,  however,  may  not  be  sufficient  to  guarantee  Q  >  0  when  fc  <  1.  The 
assumption  of  A^’s  positivity  will  play  a  crucial  role  in  finding  appropriate 
values  of  Aj,  and  C  in  numerical  examples. 

Although  an  analytical  verification  of  whether  a  particular  extremal  trajec¬ 
tory  minimizes  the  functional  is  cumbersome,  the  graph  of  this  trajectory  im¬ 
mediately  reveals  what  kind  of  an  extremal  it  is.  Indeed,  if  the  line  passing 
through  points  A  and  B  separates  the  trajectory  and  the  radar/sensor,  i.e.  the 
trajectory  moves  away  from  a  detecting  installation  (the  trajectory  is  “con¬ 
cave”),  then  it  minimizes  the  risk,  and  vise-versa,  if  the  trajectory  moves  to¬ 
wards  the  detecting  installation  (the  trajectory  is  “convex”),  then  it  maximizes 
the  risk. 

An  analytical  solution  to  equation  (11)  with  boundary  conditions  (12)  and 
constraint  (13)  is  presented  in  the  next  theorem. 


Theorem  3.3.  (analytical  solution  in  the  case  of  a  single  radar  or  sensor).  An 
analytical  solutionfor  nonlinear  first-order  differential  equation  (11)  with  con¬ 
ditions  (12),  (13)  is  given  by  the  following  quadrature 


(18) 

where  v*{p,X  l,  C)  is  a  positive  root  of  the  following  algebraic  equation  (quar- 
tic  equation) 

f(v)  =  {Cv  -  XlP^)  +  (1  -  «2)p2  -V  =  Q,  (19) 

and  unknown  constants  Xl  and  C  are  found  from  the  conditions 

(1)  /  v*(p,  Ax,,  C)  dtp  =  h  and  =  V’B  if  the  length  constraint  is 

V’A 

active; 

(2)  Ai  =  0  and  ipips)  =  i’B  if  the  length  constraint  is  inactive. 

Proof  The  main  technique  for  solving  any  first-order  differential  equation  an¬ 
alytically  is  to  explicitly  express  the  derivative  of  an  unknown  function.  By 
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introducing  an  auxiliary  function 

V = po) 

we  reduce  (11)  to  the  algebraic  equation  (19)  with  respect  to  v,  which  is  a 
particular  case  of  the  following  quartic  equation 

^2(n-2)  ^Cv  -  ALp2)2  («2^2  +  _  ^2  ^  q 

Explicit  analytical  expressions  for  four  roots  of  any  quartic  equation  may  be 
presented  by  Cardan’s  (Cardan-Ferrari’s)  formulas.  This  is  a  crucial  point  in 
obtaining  an  anal)dical  solution  for  the  differential  equation  (11).  Due  to  the 
cumbersome  form  of  the  expressions  for  the  roots  of  equation  (19)  we  do  not 
present  them  here.  Suppose  that  v*{p,  Al,  C)  is  a  root  for  (19),  then  according 
to  (20),  derivative  is  expressed 


p'^  =  ±^/{v*{p,XL,C)f-p^ 

which  leads  to  a  quadrature  expression  for  'tp  =  ip{p) 


'•Pip)  -  ^  f 


dp 


+  D. 


Excluding  constant  D  based  on  boundary  conditions  •tp{pA)  =  and 
^(pb)  =  V’B.  this  quadrature  is  reduced  to  the  form  of  (18). 

Note  a  root  for  (19)  depends  on  values  of  A^  and  C.  Which  root  should  be 
chosen  with  respect  to  Ax,  and  C  and  what  are  the  estimates  for  Xl  and  C  are 
the  subject  of  the  next  theorem.  D 

The  quadrature  ( 1 8)  is  considered  to  be  an  analytical  solution,  since  the  roots 
of  the  quartic  equation  (19)  may  be  expressed  by  Cardan’s  (Cardan-Ferrari’s) 
formulas  analytically.  There  are  two  special  cases  when  the  quadrature  (18)  is 
simplified. 


Example  1  (the  optimization  problem  without  a  constraint  on  trajectory 
length).  The  first  case  corresponds  to  the  optimization  problem  without  a  con¬ 
straint  on  trajectory  length,  in  this  case  an  optimal  trajectory  is  presented  by 
Rhodenea  (rose  function) 


p(V’)  =  C  sin 


(< 


■ii’.+  D) 


(21) 


where  D  is  a  constant  D  =  arcsin  {Cp^  —  ip  a- 
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Detail.  In  the  case  without  a  constraint  on  trajectory  length,  Ax,  =  0.  Conse¬ 
quently,  the  only  one  feasible  root  for  (19)  satisfying  t;*  >  0  is 

,  _  yi  -  (1  - 

Its  substitution  into  (18)  leads  to 

ip{p)  =  ^A±  arcsin  , 

which,  being  rewritten  as  a  function  p  =  p(^).  is  reduced  to  (21). 

In  the  case  of  n  =  2,  function  (21)  represents  an  arc  of  a  circle  passing 
through  the  origin  of  the  system  of  coordinates  and  points  A  and  B  [22].  Figure 

1.3  illustrates  behavior  of  function  p(^)  =  sin^  for  parameters 

n  =  4  and  k  =  0.5, 1.0, 2.0. 


Figure  1.3.  Function  p(V')  =  sin^  ^  . 

Note  if  n  >  2  constant  C  in  (21)  can  be  determined  only  when  |  V'b  -  V’a  I  < 
min  |7r,  otherwise  a  solution  to  (1 1- 12)  without  constraint  (13)  will  be 
unbounded. 

Example  2  (the  case  of  sphere).  The  second  case  corresponds  to  the  opti¬ 
mization  problem  when  an  aircraft  is  modeled  by  a  sphere,  in  this  case  k  =  1 
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and  an  optimal  trajectory  is  presented  by  the  explicit  quadrature 


r" 

^(Ai,r”  +  l)2-C2r2(»»-i)' 


(22) 


In  the  case  of  n  =  2,  quadrature  (22)  is  reduced  to  the  elliptic  sine  [22]. 
Detail.  In  the  case  of  k  =  1,  the  root  for  (19)  is  given  by 

Qpn-2  ’ 

which  being  substituted  into  the  quadrature  (18)  reduces  it  to  (22). 

Figure  1.4  illustrates  optimal  trajectories  for  a  “spherical”  aircraft  (k  = 
1)  for  n  =  4  with  different  constraints  on  trajectory  length,  i»,  in  trajec¬ 
tory’s  plane  determined  by  points  {xA,yA)  =  (—0.25,0.25),  ixB,yB)  = 
(1.75, 0.25)  and  radar  position  (0, 0).  Figure  1.5  shows  the  same  optimal  tra¬ 
jectories  for  a  “spherical”  aircraft  («  =  1)  for  n  =  4  with  the  same  constraints 
on  trajectory  length,  Z»,  in  3D  space  with  {xa,  yA,  ^a)  =  (—0.25, 0.15, 0.2), 
(xB,yB,ZB)  =  (1.75,0.15,0.2).  Similar  figure  forn  =  2  can  be  found  in 
[22]. 


Figure  1.4.  Optimal  trajectories  for  the  case  of  “sphere”  vnth  different  constraints  on  the 
length,  U ,  in  the  trajectories’  plane. 
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We  compared  analytical  solutions  for  optimal  trajectories  in  the  cases  with 
radars  and  sensors.  Analytical  solution  for  an  optimal  trajectory  with  a  con¬ 
straint  on  trajectory  length  in  the  case  of  a  single  sensor  is  expressed  by  the 
elliptic  sine  [22].  Figure  1.6  compares  optimal  trajectories  with  the  same  con¬ 
straint  on  trajectory  length  Z*  =  3.2  for  a  “spherical”  aircraft  in  the  cases  of 
a  single  sensor  (n  =  2)  and  single  radar  (n  =  4).  As  expected,  an  optimal 
trajectory  is  more  sensitive  to  a  radar  than  to  a  sensor  within  proximity  to  an 
installation  and  vise- versa. 


Figure  1.6.  Comparison  of  optimal  trajectories  in  trajectories’  plane  for  the  cases  with  a 
single  sensor  (n  =  2)  and  single  radar  (n  =  4)  with  the  same  constraint  on  the  length,  I,  =  3.2. 

The  next  theorem  provides  some  insights  regarding  bounds  of  unknowns  v* 
and  C  facilitating  numerical  implementation  of  (18)  and  (19). 

Theorem  3.4.  (estimates  for  \l,C  and  a  feasible  root  for  the  quartic  equa¬ 
tion).  Tjf  Ax,  >  0  then  constant  C  should  satisfy 

C  e  <  ^  ’ 

0,  min iXiPA  +  ^ ^lPb  +  ^  > 

(23) 

with  =  max  equation  (19)  has 

a  unique  root  v*  within  interval  [umin,  i^max]  where 

Vmin  =  max{p,  C7~^  (XlP^  +  min  {  1,  k~^}) , 

f  y/max{C'-2p-2(”"^)  +  0}}  . 


(24) 
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^^max  =  (kXlP^  +  p  +  C/9\/max{l,K2}  _  .  (25) 

Proof.  All  estimates  for  C  and  root  v*  are  obtained  by  analyzing  equation 
(19)  with  respect  to  feasibility  of  v*  and  the  assumption  of  A^’s  positivity.  The 
first  part  establishes  bounds  for  C  depending  on  Ax,.  Using  inequality  v  = 
+  p2  >  poTv>  p,  equation  (19)  is  reduced  to  {Cv  —  Ax,p^)  - 

V  <  0  or,  equivalently,  to  Ax,p”+^  >  -  l)  v.  Applying  v>p  again, 

the  last  inequality  is  reduced  to  A£,p”  >  Cp^~'^  - 1,  which  is  rearranged  in  the 
form 

C'<Aip  +  p-('‘-').  (26) 

Since  (26)  holds  for  all  feasible  p,  we  obtain  C  <  nun  (Ax,p  + 
Expression  Ax,p  +  is  a  convex  function  with  respect  to  p,  achieving 

its  global  minimum  at  po  =  ”•  Consequently,  if  min  {p^i,  pb}  < 

po  <  max{pA,PB}  (equivalently  Ax,  €  [A£,A]J^])  then  po  is  feasible  and 

C  <  AlPo  +  Po  j  ”  .  If  po  ^  [min  {pA,  Pb}  ,  max  {pA,  Pb}] 

(or  equivalently  Ax,  ^  [-^l  >  -^b]  )  then  we  are  not  guaranteed  that  in  a  particular 
example,  function  Ax,p  +  will  achieve  its  global  minimum  at  po,  since 

Po  may  not  be  feasible.  However,  in  this  case,  at  least  the  following  weak 

estimate  should  hold  (7  <  min  |Ax,pa  +  Pa^  Ax,Pb+Pb^”  Based 
on  w  >  0  and  the  assumption  of  Ax,  >  0,  positivity  of  C  is  obvious.  Indeed, 
rewriting  (19)  in  the  form 


we  obtain  Cv  >  XlP^,  which  concludes  that  C  >  0.  This  finalizes  the  proof 
of  formula  (23). 

The  second  part  establishes  bounds  for  v*,  i.e.  interval  [umiiD  Vmax]  con¬ 
taining  a  single  root  v*.  This  part  includes  the  following  consecutive  steps. 

■  The  first  lower  estimate  for  v*  is  obtained  by  expressing  Ax,  from  (19) 
and  satisfying  the  condition  Ax,  >  0.  That  is,  from 


we  have 


max  {C7“2p-2(n-i)  q.  ^2  _  0}. 


(27) 
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Then  we  utilize  upper  and  lower  estimates  for  -y/ +  (1  —  de¬ 
pending  on  whether  k  <  1  of  k  >  1. 

■  In  the  case  of  k  <  1,  we  have  kv  <  yj +  (1  ~  k^)p^  <  v,  which 

being  applied  to  equation  (19)  reduces  it  to  (Cv  —  XlP^)  <  1  < 

pn-2  These  inequalities  give  the  upper  and  lower  esti¬ 

mates  for  V  when  /c  <  1 

C-^  [XlP^  +  <v<C-'^  (XlP^  -H  •  (28) 

Note  if  K  =  1  then  (28)  provides  an  exact  value  for  the  root  v*  = 

C'-l(Aip2+p-(n-2)). 

■  Analogously,  in  the  case  of  «  >  1,  we  use 

KV  —  p\/ —  1  <  yj K^v"^  4-  (1  —  «2)p^  <  Kt>, 

to  reduce  equation  (19)  to 

pn-2  ^Qy  _  —  pV K^  —  1^  <v<  p"“^/ct;  (Cv  —  A^p^)  , 

where  the  left  inequality  is  then  transformed  to 

V  >  p”“^t)  ^CkV  —  X^Kf?  —  Cp\/ —  1^  . 

Consequently,  we  obtain 

v>C-^(XlP^  +  k-^P-^^-'^'^),  _ 

V  <  [XlP^  -H  k“ +  Opy/\  -  k“2J  . 


Combining  inequality  v  >  p  with  (27),  (28)  and  (29)  for  both  cases  fe  <  1 
and  K>  1,  we  obtain  (24)  and  (25). 

To  prove  that  equation  (19)  has  a  single  root  in  the  interval  [wmim  ^max].  we 
show  that  the  function 

f(v)  =  p”“^  (Cv  -  XlP^)  y^/K?V^  -t-  (1  -  «^)p^  -  V, 


is  monotonically  increasing  on  [t»min)t^max]  and  /(t^min)  <  0,  /(wmax)  >  0. 
Consider 


^/(v)  =  (Cp^  +  (1  ~  K^)p^  — 


pn-2  ^(^y  _X^p2^  f^2y 

^  y/K^V^  +  (l-K^)f^  ■ 
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The  first  term  and  the  nominator  of  the  second  term  in  the  expression  of 
are  increasing  functions  with  respect  to  v.  For  v  >  Umin  the  first 
term  is  always  nonnegative  due  to  (27)  and  the  nominator  of  the  second  term 
is  always  positive  based  on  (28)  and  (29).  Consequently,  is  positive  on 

[i^min)  Vmax].  which  means  that  f(v)  is  a  monotonically  increasing  function. 

Since  Vmin  is  the  maximum  of  three  values  (24),  we  check  the  sign  of  f{v) 
for  each  of  them. 

■  The  relation  f{p)  =  p”  (C  -  XlP  -  <  0  holds  by  virtue  of 

(26). 

■  For  Vmin  =  O},  two  cases  are  con¬ 
sidered.  If  +  K^  —  1  <  0  (when  «  <  1)  then  Omin  =  0  and 

/(f>min)  =  <  0.  If  c~^ -h  -  1  >  0  then 

Vrmn  =  f  -h  k2  _  1  and 

/ (f’min)  =  C  ^  (C'f^min  ~  ~  '£’111111  =  ~C'  XlP  ■<  0. 

■  For  Omin  =  {XlP^  +  min  {l,  based  on 


-f  (1  -  «2)p2  <  V  min{l,K}  , 

we  obtain 

fiVmin)  <  Vmin  {CVmin  “  AiP^)  min  {1,  k}  -  l) 

=  Vmin  (min  {1,  K“^}  min  {1,  «}-!)=  0. 

Thus  we  established  that  /(vmin)  <  0  for  t^min  =  max  {p,  Omim  Vmin}- 
In  the  case  of  Vmax  given  by  (24),  we  use 

y/K^V^  +  (1  -  Kp'ifP  >KV-  pi/max{l,K2}  -  1, 

to  show  that 

/(Vmax)  >  (Cvinax  “  A^p^)  (/Climax  “  P\/max  {1,  -  l)  -  Vmsx 

>  p^~‘^Vmax  (C  K  Vmax  “  Cp-^/max  {1,  _  1  -  kXl  P^)  -  Vmax  =  0. 

Consequently,  we  proved  that  /(vmin)  <  0  and  /(vmax)  >  0,  which  along 
with  the  condition  of  /(v)’s  monotonicity  on  [t;min>  Vmax]  guarantee  existence 
of  only  a  single  root  for  /  ( v)  on  [vmin  >  Vmax]  •  Q 
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Example  3  (elongated  and  compressed  ellipsoids  in  the  case  of  n  —  4). 
Coordinates  of  points  A  and  B  are  the  same  in  all  examples.  In  trajectory’s 
plane  (xA.y^)  =  (-0.25,0.25),  (»B,yB)  =  (1.75,0.25),  and  in  3D  space 
{xA,yA,ZA)  =  (-0.25,0.15,0.2),  ixB,yB,ZB)  =  (1.75,0.15,0.2).  Figures 
1.7,  1.8  and  1.9  compare  optimal  trajectories  for  sphere,  elongated  and  com¬ 
pressed  ellipsoids.  Table  1.1  presents  values  for  the  optimal  risk,  Ab,  C  and 
C’s  estimate  (23)  for  all  considered  numerical  examples  with  different  k  and 

tut . 


Figure  1.7.  Optimal  trajectories  for  sphere  (/c  =  1.0)  and  elongated  ellipsoids  (/t  —  0.5, 0.1) 
forn  =  4  and  /»  =  3.2  shown  in  trajectoiy’s  plane. 


ro6/e  1.1.  Results  of  numerical  experiments;  values  of  optimal  risk,  Ai,,  C  and  C’s  estimate. 


K 

U 

Risk 

Al 

C 

4(At/3)^^ 

1.0 

2.6 

9.792116 

4.763580322 

5.369280470 

5.658081121 

1.0 

3.2 

8.421726 

1.040107422 

1.759684438 

1.807289381 

1.0 

4.0 

7.966210 

0.300707031 

0.712062456 

0.712568693 

0.1 

3.2 

0.468371 

1.451660156 

2.282777298 

2.320693404 

0.5 

3.2 

3.980716 

1.993432500 

2.908499104 

2.943880665 

2.0 

3.2 

12.464087 

0.610351562 

1.143055224 

1.211725076 

10.0 

3.2 

14.845494 

0.109076172 

0.251462743 

0.333055390 

Analyzing  optimal  trajectories  in  Figures  1.6-1. 9  and  computational  results 
in  Table  1.1  we  conclude  the  following 

■  The  optimal  risk  is  more  sensitive  to  variation  of  the  shape  of  ellipsoid 
(parameter  k),  than  to  the  variation  of  a  trajectory’s  total  length,  i*. 
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■  Optimal  trajectories  for  different  k  (especially  for  k  >  1)  are  close  to 
each  other,  which  indicates  that  a  variation  of  ellipsoid  shape  has  no 
strong  effect  on  the  geometry  of  an  optimal  trajectory. 

m  Within  proximity  to  an  installation,  an  optimal  trajectory  is  more  sen¬ 
sitive  to  a  radar-installation  than  to  a  sensor-installation  and  in  the  area 
remote  from  the  installation  the  effect  is  opposite. 


4.  Network  Flow  Optimization  Approach 

The  calculus  of  variations  approach  reduces  the  optimization  problem  to 
the  vectorial  nonlinear  differential  equation.  Obtaining  an  analytical  solution 
to  this  equation  in  the  case  of  arbitrary  number  of  detecting  installations  is 
still  an  open  issue.  Certainly,  various  gradient-based  techniques  may  solve  the 
equation  numerically.  However,  regardless  of  efficiency  of  those  techniques 
(although,  this  issue  is  also  questionable  due  to  strong  nonlinearity  of  the  equa¬ 
tion),  most  of  them  provide  only  locally  optimal  solution.  This  section  devel¬ 
ops  a  discrete  optimization  approach  generating  globally  optimal  trajectories. 

We  propose  network  flow  (NF)  optimization  approach  to  directly  solve  the 
original  problem.  This  approach  reduces  optimal  risk  path  generation  with  a 
constraint  on  the  length  to  the  Constrained  Shortest  Path  Problem  (CSPP)  for 
a  3D  network,  which  can  efficiently  be  solved  by  NF  optimization  algorithms. 
There  are  several  advantages  of  using  NF  optimization 

■  Among  all  feasible  approximated  trajectories  in  a  considered  network,  it 
finds  a  globally  optimal  one. 

■  Its  complexity  (running  time)  depends  neither  on  a  number  of  installa¬ 
tions  in  a  network  nor  on  power  n  in  the  risk  functional  (2). 

■  It  can  readily  be  generalized  for  the  case  with  an  actual-tabulated  radar 
cross-section  (RCS)  (i.e.  when  RCS  is  not  a  smooth  function). 

However,  due  to  NP-hard  nature  of  the  CSPP,  no  polynomial  algorithm 
solves  the  CSPP  exactly.  It  means  that  in  a  worst  case,  computational  time 
for  the  CSPP  will  exponentially  depend  on  the  number  of  arcs  in  a  network. 
Consequently,  coarseness  of  the  network  should  be  specified  reasonably. 

4.1.  Network  Structure 

We  assume  an  admissible  deviation  domain  for  aircraft  trajectory  to  be  an 
undirected  graph  g  =  (Af,  A),  where  Af  =  {1, . . . ,  n}  is  the  set  consisting 
of  n  nodes  and  A  is  the  set  of  undirected  arcs.  A  trajectory  (a;(.),  y(.),  2;(.))  is 
approximated  by  a  path  V  in  the  graph  G,  where  path  V  is  defined  as  a  sequence 
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of  nodes  {jo,  ji,...,  jp)  such  that  jo  =  A,  jp  =  B  and  {jk-i,  jk)  €  ^  for  all 
k  from  1  to  p.  Let  graph  ^  be  a  3D  grid  of  nodes  (rectangular  parallelepiped) 
of  rix  X  Uy  X  riz  size  with  edges  oriented  along  coordinate  axes  x,  y,  z  and 
having  Ux,  Uy  and  riz  numbers  of  unit  segments  in  each  edge,  respectively. 
Similarly,  in  2D  case  graph  ^  is  a  2D  grid  of  nodes  (rectangle)  of  x  Uy  size 
with  edges  oriented  along  coordinate  axes  x,  y  and  having  rix  and  riy  numbers 
of  unit  segments  in  each  edge,  respectively.  Structures  of  arcs  assigned  in  Q 
in  2D  and  3D  cases  are  shown  in  Figures  1.10  and  1.11,  respectively.  A  2D 
network  with  arcs  structure,  shown  in  Figure  1.10,  contains  {nx  +  1)(%  +  1) 
nodes  and 

2(S  TXx'i^y  rix 

arcs,  where  >  1  and  Uy  >  1.  In  3D  case  the  total  number  of  nodes  and 
arcs  in  an  undirected  Q  with  arcs  structure,  as  shown  in  Figure  1.1 1,  are  {rix  + 
!)(%  +  l){nz  +  1)  and 

2(49  UxTiyriz  —  8  rixTiy  —  8  rixTiz  —  8  %nz  +  nx  +  'ny  +  riz)  > 

for  rix  >  1.  %  >  1>  and  riz  >  1-  For  instance,  the  case  ofnx  =  l,riy  =  1,  and 
riz  =  1  corresponds  to  a  single  cube  with  12  liner  arcs,  12  planar  arcs  and  4 
3D  arcs  (see  Figure  1.11).  All  these  numbers  should  be  doubled  due  to  “undi¬ 
rectness”  of  the  graph.  Thus,  the  total  number  of  arcs  just  in  a  single  cube  is 
56.  Moreover,  in  order  to  provide  sufficient  amount  of  feasible  directions  for 
a  trajectory  (i.e.  to  avoid  “naive  discretization,”  sometimes  referred  to  as  the 
digitization  bias,  [18]),  we  assign  not  only  axis  and  diagonal  but  also  so-called 
“long-diagonal”  arcs  connecting  opposite  vertexes  of  any  two  neighbor  cubes 
(see  Figure  1.1 1).  However,  in  this  case,  network  structure  becomes  very  con¬ 
dense.  For  example,  a  relatively  small  3D  undirected  network  of  40  x  40  x  40 
contains  about  69,000  nodes  and  6,200,000  arcs.  It  would  be  naive  to  assume 
that  even  a  very  efficient  NF  algorithm  is  capable  to  find  a  constrained  short¬ 
est  path  in  this  network  within  seconds  (at  least  at  the  current  moment  of  the 
technological  progress).  Obviously,  a  main  task  in  this  case  is  development  of 
efficient  network  structures  with  relatively  small  numbers  of  arcs  and  nodes, 
while  preserving  flexibility  for  trajectory  approximation,  rather  than  finding  the 
most  efficient  NF  algorithm  (although  this  is  also  a  quite  legitimate  question). 
However,  the  chapter  only  partially  addresses  the  issue  of  efficient  network 
structures,  since  this  is  a  separate  subject  for  discussion,  which  calls  for  a  sep¬ 
arate  publication.  Zabarankin  et  al.  [22]  showed  that  existing  NF  algorithms 
[7]  are  quite  efficient  in  finding  a  constrained  shortest  path  in  a  network  with 
about  100,000  arcs. 

Smoothing  procedure  and  curvature  constraint.  Number  of  arcs  can  signif¬ 
icantly  be  reduced  by  smoothing  network  structure.  Starting  from  both  jo  =  A 
and  jp  =  B  nodes  along  all  directions  outgoing  from  jo  and  jp,  we  retain  only 


3  linear  arcs  in  a  cube  6  planar  arcs  in  a  cube 


4  3D  arcs  in  a  cube 


4  planar  arcs  in  two  cubes  4  3D  arcs  in  two  cubes  2  diagonal  arcs  in  two  cubes 


Figure  1.11.  Structure  of  arcs  in  every  node  in  a  3D  network. 
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those  pairs  of  arcs,  which  produce  the  angle  not  greater  than,  for  instance, 
see  Figure  1.12. 

Let  Mi  be  the  set  of  nodes  connected  to  node  i.  For  instance,  for  2D  network 
structure  shown  in  Figure  1 . 1 0,  set  A/i  for  each  i  consists  of  1 6  nodes  producing 
4  axis,  4  diagonal  and  8  long  diagonal  arcs  with  node  i.  Let  be  a  subgraph 
starting  from  node  jo  =  A  along  arc  {jo,  i),i  ^  Mj^.  Now  TA  =  {jo,  «}•  If 
j  G  Mi/ {jo}  satisfies  condition  Ojoi  •  eij  >  where  is  the  unit  vector 
along  arc  {i,  j) ,  then  j  is  added  to  the  subgraph  7A  =  TAu  { j}.  The  next  step 
is  to  examine  new  added  nodes,  i.e.  for  all  j  €  TA/ {jo,  0  check  •  Bjk  > 
k  G  Mj.  If  node  k  satisfies  this  condition  then  TA  =  TA  u  {k}  and  so  on. 
ta  is  full  when  there  is  no  nodes  left  satisfying  the  condition.  Similarly,  we 
construct  TA  for  alH  €  Mjq.  Then  the  whole  process  is  repeated  for  TA. 


network  smoothing 


Figure  1.12.  Network  smoothing  and  curvature  constraint. 


To  avoid  sharp  turns  in  aircraft  trajectory,  we  may  use  a  curvature  constraint, 
which  also  can  be  imposed  by  aircraft  control  limitations.  Analytically,  for  any 
given  triad  of  arcs,  curvature  constraint  is  expressed  by  •  ^jkjk+i  ^ 

cos  a,  where  a,  for  instance,  may  be  In  general,  a  should  be  a  function  of 
the  length  of  the  middle  arc  {jk-i,jk),  since  a  constant  constraint  on  trajectory 
curvature  may  prevent  from  obtaining  an  optimal  solution. 

Finding  a  globally  optimal  risk  path  subject  to  length  constraint  is  the  task 
for  NF  optimization.  Network  structure  smoothing  will  be  integrated  into  a  NF 
algorithm  as  a  condition  eliminating  inadmissible  arcs  in  a  network  rather  than 
implemented  as  a  separate  procedure.  Consequently,  smoothing  condition  is 
considered  now  as  an  adjustment  for  the  NF  algorithm  rather  than  property  of 
network  structure. 
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4.2.  Approximation  Scheme 

Several  schemes  for  approximation  of  optimization  problem  (5)  are  avail¬ 
able.  We  consider  one  of  them.  Let  vector  with  components  x{jk),  y{jk) 
and  z{jk)  determine  position  of  node  jk  ■  Then  a  path  T  =  {r jj 
is  a  piece-wise  linear  curve  (broken  line)  with  vertexes  at  points  rj„A:  =  l,p. 
Any  point  on  the  arc  {jk-i,jk)  can  be  defined  by  vector  rk{t)  =  (1  —t)  rj^_^  -f 
with  t  €  [0, 1].  Thus,  length  differential  ds  and  derivative  f  for  each  arc 
are 


n  = 


^jk  ^jk-1 


k  =  l,p. 


\\^jk-^jk-i\\ 

Using  approximations  for  r,  r  and  ds,  functional  (4)  and  trajectory  length 
are  presented,  respectively 


P  ^  P 

.F(r,r)  I — rjj._j||  f  L{rk{t),ik)dt  =  (30) 

ik=i  i  k=i 

k=l 

where  \\rj^  -  J|  and  ,  rj^)  are  the  length  and  risk  index  of  the  arc 

{jk-i,jk)>  respectively.  To  derive  the  formula  for  C{rji^_^,rj^.)  we  compute 
the  risk  accumulated  along  the  arc  {jk-i,jk)  from  the  radar  located  at  qj  = 
(ci,  bi,  Ci).  Substituting  rik{t)  =  rfc(t)  -  q*  into  (2),  we  have 


>  ^jk  )  _ _ _ 

“  I+i?  l|rifcWII"+^ 

t=l  0 


_  2k 

~~  T+i? 


where 


‘iJk-l  jfc 

2 


2 

dt 


(32) 


and  </>i,  €  [0,  tt]  is  the  angle  between  vectors  rj^_^  -  qi  and  -  rj^_^ 

(see  Figure  1.13),  i.e. 
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jk-ljk  SACCOS 


/  (rjfc-1  -  qj)  ■  (rjfc  -  rjfc-i)  \ 

\lkjfc-i  ~  *iill  ~ ^jk-iWJ 


Figure  1.13  illustrates  a  3D  network  for  solving  the  risk  minimization  prob¬ 
lem.  Broken  line  AB  is  a  path  in  the  area  with  the  radar,  while  is 

the  length  of  arc  {jk-i,jk)  between  nodes  jk-i  and  jk  in  this  path.  Magnitude 
4>i,jk-\jk  angle  between  vector  rji^_i  —  qt  and  arc  {jk-i,jk)  directed 
from  node  jk-i  to  node  jk- 


Figure  1.13.  3D  networic  for  solving  the  risk  minimization  problem  -  broken  line  AB  is  a 
path  of  the  aircraft. 

Integral  (32)  can  efficiently  be  approximated  by  the  Gaussian  quadrature.  If 
/(f)  is  a  bounded  smooth  function  on  [0, 1]  then  the  Gaussian  quadrature  is 

j 

3=1 

where  hj  and  tj  €  [0, 1]  are  known  for  any  given  J.  For  instance,  Table  1.2 
presents  values  tj  and  weight  coefficients  hj  for  the  Gaussian  quadrature  for 
J=16. 

Consequently,  using  a  direct  method  of  Calculus  of  Variations,  problem  (5) 
is  approximated  by 


/ 
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Table  1.2.  Values  tj  and  weight  coefficients  hj  for  the  Gaussian  quadrature  for  J  =  16. 


2 _ V _ 

1  0.048307665687738316235 

2  0.144471961582796493485 

3  0.239287362252137074545 

4  0.331868602282127649780 

5  0.421351276130635345364 

6  0.506899908932229390024 

7  0.587715757240762329041 

8  0.663044266930215200975 

9  0.7321 821 1 8740289680387 

10  0.794483795967942406963 

11  0.849367613732569970134 

12  0.896321 1 557660521 23965 

13  0.934906075937739689171 

14  0.964762255587506430774 

15  0.98561 1511 545268335400 

16  0.997263861849481563545 


hj 

0.09654008851 4727800567 
0.09563872007927485941 9 
0.093844399080804565639 
0.091 1 7387869576388471 3 
0.08765209300440381 1 1 43 
0.08331 1 924226946755222 
0.0781 93895787070306472 
0.0723457941 08848506225 
0.065822222776361 846838 
0.0586840934785355471 45 
0.0509980592633761 76196 
0.042835898022226680657 
0.03427386291 3021 4331 03 
0.025392065309262059456 
0.01 6274394730905670605 
0.00701 861 0009470096600 


min 

P  k=l 

s- 1.  E  Ikjfc  ~ 

fe=i 

rjo  =  rA,  rjj,  =  rB- 


(33) 


If  for  all  ife  =  l,p,  Tj-fc  is  variable  (not  fixed  in  nodes  of  a  network),  then,  in 

P 

the  case  of  active  constraint  E  Ikjk  ~  II  optimality  conditions  for 

k=l 

(33)  can  be  derived  by  standard  calculus 


^ik 

where  k  =  l,(p-  1),  7  is  the  Lagrange  multiplier  for  the  constraint  in  (33)  and 
3k  =  System  (34)  may  be  solved  numerically  by  a  gradient- 

based  algorithm.  However,  in  this  case,  we  most  likely  obtain  only  locally 
optimal  solution.  Moreover,  instead  of  solving  (34)  we  could  numerically  solve 
differential  equation  (7). 
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4.3.  Reduction  to  the  Constrained  Shortest 
Path  Problem 

To  formulate  (33)  as  a  network  flow  optimization  problem,  let 


and  values  R(V)  and  1{V)  define  the  total  cost  (risk)  and  weight  (length)  ac¬ 
cumulated  along  the  path  V,  respectively, 

P  p 

fc=l  fc=l 

Thus,  each  arc  {jk-idk)  e  ^  is  associated  with  its  length  and 

nonnegative  cost  .  The  path  V  is  weight  feasible  if  the  total  weight  /(T’ ) 

is  at  most  i.e.  llv)  <  h.  Consequently,  the  CSPP  is  finding  such  a  feasible 
path  V  from  point  A  to  point  B  that  minimizes  cost  R{P) 


min 

P 

s.  t. 


p 

3k 

k=l 

p 

3k  — 

fc=i 


(35) 


The  difference  between  (33)  and  (35)  is  that  (33)  still  preserves  analytical 
properties  of  the  risk  and  length,  whereas  (35)  completely  “forgets”  about  the 
nature  of  obtained  Cjj,_j  and  values.  The  CSPP  (35)  is  closely  re¬ 

lated  to  the  Shortest  Path  Problem  with  Time  Windows  (SPPTW)  and  also  to 
the  Resource  Constrained  Shortest  Path  Problem  (RCSPP),  which  uses  a  vector 
of  weights,  or  resources,  rather  than  a  scalar.  These  problems  are  solved  in  col¬ 
umn  generation  approaches  for  Vehicle  Routing  Problems  with  Time  Windows 
(VRPTW)  and  in  long-haul  aircraft  routing  problems.  Under  the  assumption 
of  cost  and  weight  integrality,  the  CSPP  was  shown  to  be  a  NP-hard  piohltm 
[8].  It  means  that  in  a  worse  case,  the  CSPP  is  solved  in  time  exponentially 
depending  on  the  number  of  arcs.  Algorithms  for  solving  the  CSPP  are  divided 
into  three  major  categories: 


■  Label-setting  algorithms  based  on  dynamic  programming  methods 


■  Scaling  algorithms 

■  Algorithms  based  on  the  Lagrangean  relaxation  approach 

The  label  setting  algorithm  is  the  most  efficient  in  the  case  when  the  weights 
are  positive  [6].  The  Lagrangean  relaxation  algorithm  is  based  on  the  subgra¬ 
dient  optimization  [2]  and  cutting  plane  [10]  methods,  and  efficient  for  solving 
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the  Lagrangean  dual  problem  of  the  CSPP  in  the  case  of  one  resource.  Scal¬ 
ing  algorithms  use  two  fully  polynomial  approximation  schemes  for  the  CSPP 
based  on  cost  scaling  and  rounding  [11].  The  first  scheme  is  a  geometric  bi¬ 
section  search  whereas  the  second  one  iteratively  extends  paths.  We  solve  the 
CSPP  (35)  by  the  Label  Setting  Algorithm  (LSA)  with  a  preprocessing  proce¬ 
dure  [8]. 


4.4.  The  Label  Settings  Algorithm  with 
Preprocessing  Procedure 

The  Preprocessing  Procedure  and  Label  Setting  Algorithm  (LSA)  are  two 
consecutive  stages  in  finding  a  constrained  shortest  path.  The  objective  of 
preprocessing  is  to  reduce  the  original  graph  by  eliminating  all  arcs  and  nodes 
such  that  any  path  containing  them  is  infeasible  or  does  not  improve  current 
cost  upper  bound.  To  discuss  the  algorithm  in  detail,  let  us  denote  the  arc’s 
nodes  jk-\  and  jk  by  i  and  j,  respectively.  For  each  node  i,  we  consider  the 
path  obtained  by  appending  the  least  cost  path  from  the  source  node  s  to  i  to  the 
least  cost  path  from  i  to  the  sink  node  t.  If  the  total  cost  accumulated  along  the 
new  path  is  at  least  the  current  cost  upper  bound,  then  the  use  of  node  i  cannot 
improve  a  known  feasible  solution.  Hence,  node  i  and  all  arcs  incident  to  it  can 
be  deleted  from  the  graph.  If  the  total  cost  is  less  than  the  upper  bound  and  the 
path  is  feasible,  then  the  upper  bound  can  be  updated  and  the  process  continues 
with  the  improved  upper  bound.  Similar,  for  each  arc  {i,j),  we  consider  the 
path  obtained  by  appending  the  least  cost  path  from  s  to  i  to  the  least  cost  path 
from  j  to  t,  via  arc  {i,j).  If  the  total  cost  accumulated  along  the  new  path 
is  at  least  equal  to  the  current  cost  upper  bound,  then  we  can  delete  arc  {i,j) 
from  the  graph.  If  the  total  cost  is  less  than  the  upper  bound  and  the  path  is 
feasible  then  the  upper  bound  can  be  updated.  The  preprocessing  procedure  is 
presented  in  the  pseudo-code  form  below. 

Preprocessing  Algorithm  for  ttie  CSPP 

Step  0:  Let  U  =  C(n  —  1)  where  C  =  max  Cij. 

(t.i> 

Step  1 :  Find  the  minimum  cost  paths  from  source  node  s  =  A  with  arc  costs 
given  by  Cij.  Let  Q^j  be  the  least  cost  path  from  s  to  j  and  aj 
be  the  cost  of  the  path:  Oj  —  i?(QA). 

If  there  is  no  path  from  s  to  the  sink  node  t  =  B  then  stop; 
the  problem  is  infeasible. 

H  ^  then  is  the  optimal  path. 

Step  2:  Find  the  minimum  cost  paths  from  all  nodes  to  t  with  arc  costs  given 
by  Cij.  Let  be  the  least  cost  path  from  j  to  t  and  /5|  be  the 
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cost  of  the  path:  13 j  = 

Step  3:  Find  the  minimum  length  paths  from  s  to  all  nodes  with  arc  lengths 
given  by  Asy .  is  the  minimum  length  path  from  s  to  j  and  a*-  is 

the  length  of  this  path:  o:^-  =  l{Q!'sj)- 

UKQ^st)  >  i*  then  stop;  the  problem  is  infeasible. 

UKQst)  ^  ^”<1  RiQst)  <  U  then  set  U  =  iJ(Qit). 

Step  4:  Find  the  minimum  length  paths  from  all  nodes  to  t  with  the  arc  lengths 
given  by  Asy.  is  the  least  length  path  from  j  to  t  and  /3j 

is  the  length  of  this  path:  /3j  =  l{Qjt)- 

Step  5:  For  all  j  €  t}  dp 

if  o^j  +  /3j-  >  I*  then  delete  node  j  and  all  arcs  incident  to  it; 

Kaj  +  >U  then  delete  node  j  and  all  arcs  incident  to  it; 
end 

Step  6:  For  all  {i,j)  e  >4  do 

if  aj  +  Asij  +  ySj  >  Z*  then  delete  {i,  j) 
else  if  af  +  Cj j  +  /3|  >  Z7  then  delete  {i,  j) 
elseifZCQy,-')  +  Asy  +  l{Q^t)  <  I*  then  U  =  a^  +  Cij  + 

end 

Step  7:  If  during  steps  5  and  6  the  graph  changed  then  goto  Step  1, 
else  set  L  =  af  and  stop. 

End. 

The  next  stage  after  the  preprocessing  procedure  is  the  Label  Setting  Algo¬ 
rithm.  The  idea  of  the  algorithm  is  to  use  a  set  of  labels  for  each  node  and 
compare  the  labels  to  one  another.  Each  label  on  a  node  represents  a  different 
path  from  node  s  to  that  node  and  consists  of  a  pair  of  numbers  representing  the 
cost  and  weight  of  the  corresponding  path.  No  labels  having  the  same  cost  are 
stored  and  for  each  label  on  a  node,  any  other  label  on  that  node  with  a  lower 
cost  must  have  a  greater  weight.  Let  Jj  be  the  index  set  of  labels  on  node  i  and 
for  each  k  e  li  let  denote  a  path  from  s  to  i  with  weight  Wj^  and  cost  Cf . 
Pair  Cf )  is  the  label  of  node  i  and  is  the  path  corresponding  to  it.  For 

two  labels  and  C^),  corresponding  to  two  different  paths  Vi 

and  Vi,  respectively,  (Wf ,  C*- )  dominates  {Wi,  Ci)  if  Wi  <  Wf ,  Cf  <  Cf , 
and  the  labels  are  not  equal.  Label  (W/',  C^)  is  efficient  if  it  is  not  dominated 
by  any  other  label  at  node  i,  i.  e.  if  {1{V),  R{V))  does  not  dominate  {Wi ,  C^) 
for  all  paths  V  from  s  to  t.  A  path  is  efficient  if  the  label  it  corresponds  to  is 
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efficient.  The  LSA  finds  all  efficient  labels  in  evety  node.  Starting  without  any 
labels  on  any  node,  except  for  label  (0,  0)  on  node  s,  the  algorithm  extends  the 
set  of  all  labels  by  treating  an  existing  label  on  a  node,  that  is,  by  extending 
the  corresponding  path  along  all  outgoing  arcs.  Let  Li  be  the  set  of  labels  on 
node  i  and  let  Tj  C  /j  index  the  labels  on  node  i,  which  have  been  treated.  The 
algorithm  proceeds  until  all  labels  have  been  treated,  i.e.  until  Ii\Ti  =  0  for 
all  i  6  F\{0- 


The  Label  Setting  Algorithm  (LSA)  with  smoothing  condition 


Step  0;  Initialization 

Run  Preprocessing  Algorithm  for  the  CSPP  to  find  U,  y0|,  /?]• 
andQ|e  VjeV\{t}. 

Set  Ls  =  {(0,  0)}  and  Li  =  iti  for  all  i  €  V\{s}. 

Initialize  L  accordingly  for  each  i  €  V. 

Set  Ti  =  ^  for  each  i  6  V. 

Step  1 :  Selection  of  the  label  to  be  treated 

If  U  {Ii\Ti)  =  0  then  stop;  all  efficient  labels  have  been  generated, 
iev 

Else  choose  i  eV  and  k  E  Ii\Ti  so  that  Wf  is  minimal. 

Step  2:  Treatment  of  label  (Wf ,  C^) 

For  all  {i,j)  €  v4dp 

ff  (e(<-i)i  •  ^ij  >  e)  I* smoothing  condition:  e  =  (♦—  i)  is 

a  predecessor  node  */ 

If(Wf +  Asij+;05  <L) 

U(C^  +  Cii+Pl<U) 

K  {^i  +  Asy ,  Cf  +  Cij)  is  not  dominated 
hy{W],q)  'iqElj 

then  set  Lj  =  Lj  U  {(W/'  +  Asy,  Cf  +  Cij)} 
and  update  Ij 

If  (Wf  +  Asy  +  Z(Q^t)  <  Z,)  then  t/  =  Cf  +  Cij  + 

end 

Step  3:  Set  Tj  =  Ti  U  {A:},  goto  to  Step  1. 


End. 
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5.  Numerical  Experiments 

Zabarankin  et  al.  [22]  demonstrated  efficiency  of  the  NF  optimization  ap¬ 
proach  in  optimal  trajectory  generation  in  2D  space  with  arbitrary  number  of 
sensors  (n  =  2)  for  the  case  of  sphere  (k  =  1).  This  section  tests  complexity  of 
the  CSPP  (35)  and  flexibility  of  the  proposed  3D  network  structure  in  optimal 
trajectory  generation  with  variable  aircraft  RCS  in  3D  space.  We  computed 
discrete  solutions  using  i2D  and  3D  networks  for  the  same  data  considered  in 
Examples  2  and  3  in  the  case  of  a  single  radar  (n  =  4).  Radar  position,  co¬ 
ordinates  of  points  A  and  B  and  testing  values  for  a  constraint  on  trajectory 
length  are  exactly  the  same.  In  3D  case,  the  CSPP  was  solved  by  the  LSA  with 
and  without  smoothing  condition.  All  calculations  were  conducted  using  a  PC 
with  Xeon  3.08  GHz  and  3.37  Gb  of  RAM. 

We  used  2D  and  3D  networks,  with  structures  as  shown  in  Figures  1.10  and 
1.11,  to  compare  discrete  optimization  trajectories  with  analytical  ones  in  the 
case  of  a  single  radar.  2D  network  is  a  special  case  of  3D  one  with  =  0.  It 
tests  discrete  solutions  in  the  trajectories’  plane  determined  by  points  A,  B  and 
(0, 0, 0)  (radar  position).  We  associate  nodes  of  a  3D  graph  with  integer  vectors 
(i,  j,  k)  fomiing  3D  integer  grid.  Consequently,  the  set  of  arcs  lengths  in  an 
integer  grid  with  structure  as  shown  in  Figure  1.11  is  {l,  V2,  \/3,  ^/E,  VE,  3}. 
In  assigning  real  arc  length,  all  these  values  are  scaled  by  an  appropriate  co¬ 
efficient  depending  on  actual  size  of  a  network.  In  order  to  reduce  compu¬ 
tational  time,  we  approximated  arcs  lengths  by  the  set  of  integer  numbers 
{1000, 1414, 1732, 2236, 2449, 3000}.  In  this  case,  the  scaling  coefficient  is 
adjusted  correspondingly.  Finding  a  constrained  shortest  path  in  a  network 
with  the  integer  lengths  of  arcs  is  approximately  1.5  times  faster  than  the  same 
procedure  with  real  lengths  of  arcs.  However,  due  to  the  integer  approximation 
of  arcs  lengths,  the  actual  length  of  a  constrained  shortest  paths  in  a  network 
may  be  slightly  greater  than  assigned  length  constraints,  while  the  correspond¬ 
ing  optimal  risk  value  may  be  lesser  than  the  one  obtained  by  analytical  so¬ 
lutions  approach.  In  tables  presenting  results  of  network  flow  optimization, 
optimal  risk  values,  obtained  by  network  optimization  and  inconsistent  with 
“true”  ones  in  the  discussed  sense,  are  marked  by  ^  symbol. 

We  calculated  constrained  shortest  paths  depending  on  ellipsoid  shape  (pa¬ 
rameter  k)  and  length  constraint,  and  compared  optimal  risk  values  R2D  and 
Rsd  (in  2D  and  3D  network  optimization,  respectively)  with  the  “true”  ones, 
obtained  by  analytical  solutions  approach.  Values  of  I2D  and  Izd  are  lengths 
of  constrained  shortest  paths  in  2D  and  3D  cases,  respectively.  We  were  in¬ 
terested  in  the  following  parameters:  number  of  nodes  left  after  preprocessing, 
Nprep,  cost  Upper  bound  in  preprocessing,  U,  preprocessing  time,  Tprep,  num¬ 
ber  of  labels  treated  in  the  LSA,  Niabeis^  and  running  time  of  the  LSA,  Tlsa, 
measured  in  seconds.  All  these  parameters  are  helpful  in  evaluating  perfor- 
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mance  of  the  discrete  optimization  approach.  We  analyzed  the  impact  of  using 
smoothing  condition  on  accuracy  of  discrete  solutions  and  LSA  running  time 
in  3D  case. 


5.1.  2D  Network  Optimization  in  the  Case  of  a 
Single  Radar 

To  calculate  2D  optimal  trajectories  for  the  cases  considered  in  Examples  2 
and  3,  we  used  a  2D  squared  graph  with  the  following  parameters 


Size  of  the  graph 

Tlx  X  Tly 

Length  of  axis  arcs 
Number  of  nodes 
Number  of  arcs 
Radar  position 
Point  A 
Point  B 


2.3  X  2.3 
46  X  46 
2.3/46  =  0.05 
(46  +  1)^  =  2209 
33672 
(0,0) 

(-0.25,0.25) 

(1.75,0.25) 


Numerical  results  of  2D  network  optimization  for  different  values  of  k  and 
are  presented  in  Tables  1.3  and  1.4.  Figures  1.14  and  1.15  compare  analytical 
and  discrete  optimization  solutions  in  2D  space  with  parameters;  a)  n  =  4,  k  = 
1.0,  h  =  2.6,  3.2, 4.0;  and  b)  n  =  4,  «  =  0.1, 2.0,  U  =  3.2,  respectively.  The 
smooth  curves  are  the  optimal  trajectories  obtained  by  the  analytical  approach 
and  the  nonsmooth  curves  are  those  obtained  by  solving  the  CSPP.  Analytical 
and  corresponding  discrete  optimization  trajectories  are  close  to  each  other, 
which  validates  both  approaches.  Note  for  the  case  of  k  =  10.0,  values  of 
optimal  risk  for  disctete  trajectories  in  Tables  1.3,  1.4,  1.5,  1.6  and  1.7  are 
lesser  than  the  risk  value  for  the  corresponding  analytical  solution.  Due  to 
integer  approximation  of  arcs  lengths,  the  total  length  of  those  paths  are  greater 
than  their  integer  representations. 


Table  1.3.  Results  of  2D  network  preprocessing:  single  radar. 


K 

u 

True  Risk 

u 

JV 

Tpr,  sec 

1.0 

2.6 

9.792116 

91.336509 

837 

0.359 

1.0 

3.2 

8.421726 

91.336509 

1260 

0.313 

1.0 

4.0 

7.966210 

35.693792 

1703 

0.500 

0.1 

3.2 

0.468371 

15.942982 

1281 

0.484 

0.5 

3.2 

3.980716 

66.836609 

1265 

0.329 

2.0 

3.2 

12.464087 

12.858975 

764 

0.547 

10.0 

3.2 

14.845494 

14.84242^ 

42 

0.359 
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Table  1.4.  Results  of  2D  network  optimization  with  LSA:  single  radar. 


K 

i. 

True  Risk 

hn 

R2D 

l^labels 

TlsAj  sec 

1.0 

2.6 

9.792116 

2.592 

9.983805 

85461 

0.609 

1.0 

3.2 

8.421726 

3.1995 

8.504245 

224524 

2.047 

1.0 

4.0 

7.966210 

3.9933 

8.004489 

423056 

4.750 

0.1 

3.2 

0.468371 

3.199 

0.488162 

329517 

3.641 

0.5 

3.2 

3.980716 

3.1932 

4.063807 

292594 

2.922 

2.0 

3.2 

12.464087 

3.1958 

12.518963 

60619 

0.406 

10.0 

3.2 

14.845494 

3.1958 

14.84242+ 

0 

0 

Figure  I.I4.  Comparison  of  analytical  and  discrete  optimization  trajectories  for  the  case  of 
sphere  (/c  =  1.0),  n  =  4  and  different  length  constraints,  U,  in  trajectories’  plane. 
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Figure  1.15.  Comparison  of  analytical  and  discrete  optimization  trajectories  for  elongated 
{k  —  0.1)  and  compressed  (k  =  2.0)  ellipsoids  for  n  =  4  and  the  same  constraint  on  the  length 
U  =3.2  in  trajectories’  plane. 


5.2.  3D  Network  Optimization  in  the  Case  of  a 
Single  Radar 

For  3D  network  optimization  with  the  data  from  Examples  2  and  3,  we  used 
a  parallelepiped  graph  with  parameters 


Size  of  the  graph 

Tlx  Xriy  X  Tly 

Length  of  axis  arcs 
Number  of  nodes 
Number  of  arcs 
Radar  position 
Point  A 
Point  B 


2.3  X  1.0  X  1.25 
46  X  20  X  25 

2.3/46  =  1.0/20  =  1.25/25  =  0.05 
(46  +  1)(20  -f  1)(25  +  1)  =  25662 
2213062 
(0,0,0) 

(-0.25,0.15,0.2) 

(1.75,0.15,0.2) 


Numerical  results  of  3D  network  optimization  with  and  without  network 
structure  smoothing  for  different  values  of  k  and  Z*  are  presented  in  Tables  1.5, 
1.6  and  1.7.  Figures  1.16,  1.17  and  1.18  compare  the  analytical  and  discrete 
optimization  solutions  in  3D  space  for  the  following  sets  of  parameters;  a)  n  = 
4,  K  =  1.0,  /*  =  2.6,  3.2,  4.0;  b)  n  =  4,  k  =  0.1,  h  =  3.2;  and  c)  n  =  4, 
K  =  2.0,  Z*  =  3.2,  respectively. 
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Table  1.5.  Results  of  3D  network  preprocessing:  single  radar. 


K 

i. 

True  Risk 

U 

Npr 

Tpx'f  sec 

1.0 

2.6 

9.792116 

91.336509 

11518 

9.516 

1.0 

3.2 

8.421726 

91.336509 

22521 

14.234 

1.0 

4.0 

7.966210 

8.902138 

21175 

18.031 

0.1 

3.2 

0.468371 

15.942982 

22598 

16.094 

0.5 

3.2 

3.980716 

66.836609 

22553 

14.875 

2.0 

3.2 

12.464087 

12.969952 

15543 

22.187 

10.0 

3.2 

14.845494 

14.840356^ 

2873 

12.781 

Table  1.6.  Results  of  3D  network  optimization  with  LSA:  single  radar. 


K 

u 

True  Risk 

ho 

Rzd 

^labels 

TLSAjSec 

1.0 

2.6 

9.792116 

2.5998 

10.251423 

1936613 

766.704 

1.0 

3.2 

8.421726 

3.19805 

8.525182 

6644608 

6597.750 

1.0 

4.0 

7.966210 

3.9997 

8.040955 

8066613 

7252.609 

0.1 

3.2 

0.468371 

3.1987 

0.554746 

9930869 

11519.922 

0.5 

3.2 

3.980716 

3.1963 

4.11088 

8427875 

9069.281 

2.0 

3.2 

12.464087 

3.19895 

12.529767 

2529266 

1220.188 

10.0 

3.2 

14.845494 

3.1953 

14.83876^ 

132135 

8.031 

Table  1.7.  Results  of  3D  network  optimization  with  LSA  &  smoothing  condition:  single 
radar. 


« 

i. 

True  Risk 

IZD 

Rzd 

^labels 

Tlsa,  sec 

1.0 

2.6 

9.792116 

2.5998 

10.251423 

1080091 

299.735 

1.0 

3.2 

8.421726 

3.19645 

8.526340 

4841651 

4286.844 

1.0 

4.0 

7.966210 

3.9981 

8.041533 

6849799 

5937.016 

0.1 

3.2 

0.468371 

3.19925 

0.55662 

5262314 

5066.953 

0.5 

3.2 

3.980716 

3.1986 

4.123165 

5091681 

4715.875 

2.0 

3.2 

12.464087 

3.1998 

12.530477 

2015458 

846.672 

10.0 

3.2 

14.845494 

3.1953 

14.83876+ 

76150 

4.516 

Radar 

-0.5 

Figure  1.16.  Comparison  of  analytical  and  discrete  optimization  trajectories  for  sphere  («  = 
1.0),  n  =  4  with  different  length  constraints,  /*,  in  3D  space. 


Figure  1.17.  Comparison  of  analytical  and  discrete  optimization  trajectories  for  elongated 
ellipsoid  «  =  0.1  and  parameters  n  =  4,  /*  =  3.2  in  3D  space. 
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Figure  1.18.  Comparison  of  analytical  and  discrete  optimization  trajectories  for  compressed 
ellipsoid  k  =  2.0  and  parameters  n  =  4,  /♦  =  3.2  in  3D  space. 


5.3.  3D  Network  Optimization  in  Cases  with 
Two  and  Three  Radars 

This  section  analyzes  impact  of  variable  RCS  in  the  case  of  several  radars 
on:  a)  geomety  of  optimal  trajectories;  and  b)  performance  of  the  discrete 
optimization  approach.  For  optimal  trajectory  generation  in  the  cases  with  two 
and  three  radars,  we  used  the  same  3D  network  of  2.3  x  1.0  x  1.25  with  the 
corresponding  integer  grid  nx  x  %  x  ny  =  46  x  20  x  25  and  the  following 
data  for  radars  positions  and  staring  and  ending  trajectory  points 


Radar  1  = 

(1,0,0) 

Point  A 

=  (0, 0.5,0) 

Radar  2  = 

(0.5, 1,0) 

Point  B 

=  (2, 0.5,0) 

Radar  1  = 

(1,0,0) 

Point  A 

=  (0,0.75,0) 

Radar  2  = 

(0.5,1.25,0) 

Point  B 

=  (2,0.75,0) 

Radar  3  = 

(1.5, 1,0) 

Numerical  results  of  3D  network  preprocessing  and  optimization  in  the 
cases  of  two  and  three  radars  with  and  without  network  structure  smoothing 
for  different  values  of  k  and  the  same  constraint  on  the  length,  Z*  =  3.2,  are 
presented  in  Tables  1.8,  1.9,  1.10  and  1.11.  Figures  1.19-1.24  illustrate  dis¬ 
crete  optimization  trajectories  in  3D  space  with  two  and  three  radars  for  the 
following  parameters:  n  =  4,  Z*  =  3.2,  k  =  0.1, 1.0, 2.0. 
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Table  1.8.  Results  of  3D  network  preprocessing:  two  radars 


K 

U 

Npr 

Tpr,  sec 

0.1 

3.2 

3.086443 

20342 

16.547 

1.0 

3.2 

15.367339 

19267 

12.766 

2.0 

3.2 

7.535904 

15968 

16.828 

Table  1,9.  Results  of  3D  network  optimization:  two  radars  (*  smoothing  is  used) 


K 

i. 

IzD 

Rsd 

Nlabels 

Tlsa,  sec 

0.1 

3.2 

3.19885 

0.921916 

9993019 

9837.906 

1.0 

3.2 

3.1993 

4.891124 

6737166 

5392.109 

2.0 

3.2 

3.1993 

4.320838 

3617699 

2010.750 

*0.1 

3.2 

3.1993 

0.922975 

6428635 

5556.891 

*1.0 

3.2 

3.1993 

4.891411 

5095577 

3708.437 

*2.0 

3.2 

3.19785 

4.330046 

2818012 

1335.859 

Table  1.10.  Results  of  3D  network  preprocessing:  three  radars 


K 

i. 

U 

Npr 

Tpr,  sec 

0.1 

3.2 

19.374493 

20766 

14.718 

1.0 

3.2 

28.740118 

19660 

16.641 

2.0 

3,2 

22.262851 

18377 

15.625 

Table  1.11.  Results  of  3D  network  optimization:  three  radars  C^smoothing  is  used) 


K 

1. 

IzD 

Rzd 

Nlahels 

Tlsa,  sec 

0.1 

3.2 

3.19995 

1.635309 

7898013 

6733.031 

1.0 

3.2 

3.19905 

9.073111 

5623460 

3998.969 

2.0 

3.2 

3.19905 

8.62298 

3982280 

2285.640 

*0.1 

3.2 

3.19995 

1.639800 

5296123 

4056.406 

•1.0 

3.2 

3.197 

9.094264 

4413319 

2916.281 

*2.0 

3.2 

3.19905 

8.62298 

3155044 

1634.031 
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Figure  1.19.  Optimal  trajectories  in  the  case  of  two  radars  for  compressed  ellipsoid  (k  — 
2.0),  sphere  (k  =  1.0)  and  elongated  ellipsoid  (k  =  0.1)  with  the  same  length  constraint, 
U  =  3.2. 


Figure  1.20.  Front  view:  optimal  trajectories  in  the  case  of  two  radars  for  compressed  el- 
lipsoid  (k  =  2.0),  sphere  (k  =  1.0)  and  elongated  ellipsoid  (k  =  0.1)  with  the  same  length 
constraint,  U  =  3.2. 
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Figure  1.21.  View  from  above:  optimal  trajectories  in  the  case  of  two  radars  for  compressed 
ellipsoid  (k  =  2.0),  sphere  (k  =  1.0)  and  elongated  ellipsoid  (k  =  0.1)  with  the  same  length 
constraint,!*  =  3.2. 


Figure  1.22.  Optimal  trajectories  in  the  case  of  three  radars  for  compressed  ellipsoid  (k  = 
2.0),  sphere  (k  =  1.0)  and  elongated  ellipsoid  («  =  0.1)  with  the  same  length  constraint, 
L  =  3.2. 
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Figure  1,23.  Front  view:  optimal  trajectories  in  the  case  of  three  radars  for  compressed  el¬ 
lipsoid  (k  =  2,0),  sphere  (k  =  1.0)  and  elongated  ellipsoid  (k  —  0.1)  with  the  same  length 
constraint,  U  =  3.2. 


Figure  1,24,  Side  view:  optimal  trajectories  in  the  case  of  three  radars  for  compressed  el¬ 
lipsoid  (k  =  2.0),  sphere  («  =  1.0)  and  elongated  ellipsoid  («  =  0.1)  with  the  same  length 
constraint,  U  =3.2. 
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■  There  is  no  strong  correlation  between  LSA  running  time  and  number  of 
radars;  depending  on  ellipsoid  shape  it  may  decrease  (k  =  1.0),  increase 
(k  =  2.0)  or  variate  (k  =  1.0). 

■  Running  time  of  preprocessing  procedure  is  always  small  (10-20  sec 
in  3D  case),  which  in  most  testing  examples  is  less  than  2%  of  total 
computational  time,  and  indicates  no  predictive  power  for  LSA  running 
time.  However,  number  of  nodes  leA  after  preprocessing  is  helpful  in 
evaluating  expected  LSA  running  time.  Also,  testing  examples  suggest 
that  LSA  running  time  may  linearly  depend  on  number  of  treated  labels. 
Although  this  number  is  known  only  after  the  algorithm  stops,  it  can  be 
used  as  a  reference  value  for  another  run. 

■  In  both  examples  with  several  radars,  optimal  trajectories  for  different 
values  of  parameter  k  (ellipsoid  shape),  subject  to  the  same  constraint 
on  trajectory  length,  are  again  close  to  each  other  (the  same  phenomenon 
was  observed  by  comparing  analytical  trajectories  in  the  case  of  a  single 
radar). 


Tlsa>  sec 


Figure  1.25.  Dependence  of  LSA  running  time  on  the  shape  of  ellipsoid,  k  (3D  network, 
single  radar):  curve  “1”  -  no  smoothing,  curve  “2”  -  smoothing  is  used. 

Figure  1.25  shows  dependence  of  LSA  running  time  on  the  shape  of  ellip¬ 
soid,  K,  in  the  case  of  3D  space  and  a  single  radar  with  and  without  network 
smoothing.  The  excessive  running  time  for  the  LSA  in  the  case  of  very  elon¬ 
gated  ellipsoids,  (k  «  1,  ellipsoid  with  k  =  0.1  is  almost  a  needle),  can  be 
explained  by  lowest  risk  accumulations  in  directions  radial  to  a  radar,  which 
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are  ones  producing  greatest  total  lengths  from  point  A  to  point  B.  This  com¬ 
plicates  comparison  of  labels  in  risk  minimization,  while  balancing  length  con¬ 
straint.  This  idea  is  supported  by  the  fact  that  the  running  time  in  generating 
optimal  trajectory  for  a  compressed  ellipsoid  with  k  =  10.0  is  just  several  sec¬ 
onds.  In  this  case,  because  of  compressed  geometry  (for  instance,  a  disk  flying 
along  its  axis  of  symmetry),  the  risk  of  detection  accumulates  slower  in  direc¬ 
tions  transversal  to  a  radar  and  those  direction  are  the  ones  producing  lowest 
total  lengths  from  point  A  to  point  B.  It  is  worth  to  mention  that  for  small 
values  of  k,  network  smoothing  reduces  LSA  running  time  more  efficiently. 

Running  time  analysis.  Figure  1.26  illustrates  dependence  of  LSA  mnning 
time,  TiSA.  oo  number  of  labels  treated,  Niateist  in  a  3D  network  for  a  sin¬ 
gle  radar  and  various  k  and  Z,  with  and  without  smoothing.  Variations  of  k 
and  Z*  have  no  strong  effect  on  LSA  running  time  because  of  the  preprocess¬ 
ing  procedure.  The  running  time  is  almost  linearly  depends  on  the  number  of 
labels  treated,  which,  in  turn,  is  a  function  of  number  of  nodes  left  after  pre¬ 
processing,  Npr,  smoothing  condition  and  cost  upper  bound,  U ,  obtained  in 
preprocessing.  Figure  1.27  shows  strong  correlation  between  ATjafteis  aiid  Npr 
plotted  for  all  k  and  Z*  in  the  case  of  a  single  radar.  While  Niabeis  is  uniformly 
reduced  by  smoothing  condition,  it  may  be  quite  different  for  the  same  value 
of  Npr  because  of  different  cost  upper  bounds  obtained  in  preprocessing.  Ob¬ 
viously,  the  closer  U  to  optimal  cost  is,  the  lower  number  of  labels  will  be 
treated.  According  to  results  presented  in  Tables  1.3,  1.5,  1.9  and  1.11,  cost 
upper  bounds  are  not  close  enough  to  “true”  risk  values.  This  fact  suggests  to 
develop  preprocessing  procedures  obtaining  more  accurate  cost  upper  bounds. 
Such  preprocessing  may  be  based  on  Lagrange  relaxation  [10, 1 1]. 


7.  Conclusions 

We  developed  a  three-dimensional  deterministic  model  for  routing  ah  air¬ 
craft  with  a  variable  radar  cross-section  (RCS)  in  a  threat  environment.  The 
threat  is  associated  with  the  risk  of  detection  by  radars,  sensors  or  surface  air 
missiles.  To  investigate  dependence  of  the  risk  of  detection  on  variable  RCS, 
we  model  the  aircraft  by  a  symmetrical  ellipsoid  with  the  axis  of  symmetry  ori¬ 
enting  trajectory  direction.  The  model  considers  the  risk  of  detection  to  be  the 
sum  of  risks  from  all  installations  in  the  area  of  interest,  where  the  risk  to  be 
detected  by  a  particular  installation  is  proportional  to  the  area  of  ellipsoid  pro¬ 
jection  and  reciprocal  to  the  -power  of  the  distance  between  the  aircraft  and 
this  particular  installation.  We  developed  analytical  and  discrete  optimization 
approaches  for  solving  the  risk  minimization  problem  subject  to  a  constraint 
on  trajectory  length. 

The  analytical  approach,  based  on  calculus  of  variations,  reduces  the  origi¬ 
nal  problem  to  solving  the  vectorial  nonlinear  differential  equation.  We  derived 
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Figure  1.86.  LSA  running  time  versus  number  of  labels  treated:  3D  network,  single  radar. 
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Figure  1.27.  Number  of  labels  treated  versus  number  of  nodes  left  after  preprocessing  (3D 
network,  single  radar):  curves  “1”  and  “2”  correspond  to  LSA  and  LSA  with  smoothing,  re¬ 
spectively. 
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this  equation  based  on  a  general  form  of  the  risk  functional  in  the  case  with  an 
arbitrary  number  of  passive  or  active  installations.  For  the  case  of  a  single 
installation,  arbitrary  ellipsoid  shape  and  any  n,  we  obtained  an  analytical  so¬ 
lution  to  the  vectorial  differential  equation,  which  is  expressed  by  a  quadrature. 
According  to  numerical  experiments  based  on  the  analytical  solutions,  we  con¬ 
clude  that 

■  Complexity  of  solving  the  vectorial  differential  equation  analytically  is 
primarily  determined  by  the  number  of  installations  in  the  area  of  interest 
and  is  not  affected  by  the  type  of  an  installation  (radar  or  sensor). 

■  In  the  case  of  a  single  installation: 

-  An  analytical  solution  is  expressed  by  a  quadrature  and  correspond¬ 
ing  optimal  trajectory  is  a  planar  curve  within  the  plane  defined  by 
starting  and  hnishing  trajectory’s  points  and  the  radar’s  position. 

-  The  model  with  constant  RCS  (’’spherical”  aircraft,  k  =  1.0)  es¬ 
sentially  simplifies  obtaining  an  analytical  solution  and  its  further 
numerical  analysis. 

m  The  optimal  risk  is  more  sensitive  to  the  variation  of  ellipsoid  shape  than 
to  the  variation  of  trajectory  total  length. 

■  Optimal  trajectories  for  different  k  (especially  for  k  >  1)  are  close  to 
each  other,  which  indicates  that  a  variation  of  ellipsoid  shape  has  no 
strong  effect  on  the  geometry  of  an  optimal  trajectory. 

m  Within  proximity  to  an  installation,  an  optimal  trajectory  is  more  sensi¬ 
tive  to  a  radar-installation  than  to  a  sensor-installation  and  in  area  remote 
from  the  installation  the  effect  is  opposite. 

Obtaining  an  analytical  solution  to  the  vectorial  differential  equation  in  the 
case  with  an  arbitrary  number  of  installations  is  still  an  open  issue.  However, 
availability  of  an  analytical  solution  in  the  case  of  a  single  installation  signifi¬ 
cantly  facilitates  conceptual  understanding  the  impact  of  variable  RCS  on  the 
geometry  of  optimal  trajectories  and  testing  discrete  optimization  approaches. 

To  address  optimal  trajectory  generation  in  3D  space  in  the  case  of  vari¬ 
able  RCS  and  arbitrary  number  of  radars,  we  developed  discrete  optimization 
approach  based  on  network  flow  optimization.  Approximating  the  area  of  in¬ 
terest  by  a  3D  network  with  a  flexible  structure  and  presenting  a  trajectory  by 
a  path  in  this  network,  NF  optimization  reduced  optimal  risk  path  generation 
with  a  constraint  on  trajectory  length  to  the  Constrained  Shortest  Path  Problem 
(CSPP).  We  suggested  to  solve  the  CSPP  by  Label  Setting  Algorithm  (LSA) 
with  network  smoothing,  which  is  considered  as  an  adjustment  for  the  algo¬ 
rithm  rather  than  a  property  of  network  structure.  This  condition,  intended  for 
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preserving  trajectory  smoothness  and,  as  a  result,  eliminating  inadmissible  arcs 
in  the  network,  can  be  used  as  a  necessary  constraint  in  trajectory  generation. 
We  tested  NF  optimization  approach  for  2D  and  3D  networks  with  and  with¬ 
out  smoothing  condition,  wi^  various  ellipsoid  shapes,  several  constraints  on 
trajectory  length  in  the  cases  with  one,  two  and  three  radars.  Based  on  compu¬ 
tational  results  of  testing  examples,  we  made  the  following  conclusions 

m  In  the  case  of  a  single  radar,  all  optimal  trajectories  obtained  by  the  dis¬ 
crete  approach  for  various  k  and  are  sufficiently  close  to  the  corre¬ 
sponding  analytical  trajectories. 

■  Network  smoothing  condition  reduces  LSA  running  time  by  1. 5-2.5 
times,  while  preserving  accuracy  of  optimal  trajectories. 

m  LSA  generates  2D  discrete  trajectories  within  a  few  seconds. 

■  LSA  running  time  is  extremely  sensitive  to  the  shape  of  ellipsoid;  in  3D 
case,  it  varies  from  5  to  5000 sec  for  0.1  <  k  <  10.0. 

■  In  testing  examples  with  two  and  three  radars,  optimal  trajectories  with 
the  same  constraint  on  the  length  but  different  ellipsoid  shapes  are  rela¬ 
tively  close  to  each  other,  which  suggests  that  in  general,  ellipsoid  shape 
has  no  strong  effect  on  the  geometry  of  an  optimal  trajectory. 

■  Running  time  of  the  algorithm  strongly  depends  on  the  value  of  trajec¬ 
tory  length  constraint. 

This  chapter  introduced  3D  analytical  model  addressing  optimal  trajectory 
generation  with  variable  RCS  subject  to  a  constraint  on  trajectory  length.  De¬ 
veloped  anal5rtical  and  discrete  optimization  approaches  are  just  a  first  step  in 
solving  the  proposed  model  rather  than  exhaustive  answer  to  this  matter.  In  3D 
case,  other  NF  algorithms  for  solving  the  CSPP  as  well  as  other  approximation 
schemes  for  the  original  risk  minimization  problem  may  be  addressed. 


Appendix:  Minimization  of  a  Functional  with  Non- 
holonomic  Constraint  and  Movable  End  Point 

This  section  reduces  necessary  conditions  for  minimization  of  a  functional  with  a  nonholo- 
nomic  constraint  and  a  moveable  end  point  to  a  vectorial  nonlinear  differential  equation.  This 
equation  plays  a  central  role  in  solving  (5)  in  the  case  of  a  single  radar.  We  consider  the  follow¬ 
ing  general  formulation 

min$(r,f,/),  (l.A.l) 

r 

I 

$(r,r,i)=  L(r(s),r(s))ds, 


0 


(1.A.2) 
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(1.A.3) 

(1.A.4) 

(1.A.5) 


r(0)  =  ri,  r(/)  =  r2, 

(/?(r(s))  =  0, 

1<U, 

where  r(5)  =  (x(s),  y{s),  z{s)),  r{s)  —  {x{s),y{s),  z{s)). 

A  necessary  condition  for  the  existence  of  a  functional  extremum  requires  the  total  variation 
of  the  functional  to  be  equal  to  zero.  However,  because  of  constraint  (1.A.5),  (1.A.1H1-A.5) 
is  the  problem  with  the  movable  end  point,  r(/),  which  means  that  variation  of  the  total  curve 
length,  Z,  is  not  zero.  Note  variations  6x,  Sy  and  ^2;  are  dependant  by  virtue  of  nonholonomic 
constraint  (1.A.4).  Lagrange  multiplier  method  is  used  to  separate  differential  expressions  in 
the  functional  variation. 

Let  6r  =  (Sx^Sy^Sz),  6t  =  {dx,Sy,  i),  ^  =  (||,  fj)  and  ^  •^)* 

Applying  the  Lagrange  multiplier  method  to  problem  (1.A.1)-(1.A.2)  with  constraint  (1.A.4), 

the  functional  (1  .A.2)  is  rewritten  as  #(L,  0,  A,  Z)  =  (L(r,  r)  +  A(5)  <^(r))  ds.  By  definition 

0 

the  variation  of  this  functional  is 


=  {SL{r,  r)  +  A  ^<^(r)  +  ^ A)  cZs  +  (L  +  A  0)  |  SI 
0 

=  ^  -h  •  <Jr  +  A  If  •  (5r  +  <^(5A  ds -{- {L  + X(l>)\g^iSl 

0 

=  I  .5v  +  <l>5Xds+  t  +  Af  .5r 

+  4"  A^)|^_^  JZ. 

Note  y  0,  since  Z  is  varied  and  5  =  Z  is  not  anymore  a  boundary  point.  Based  on 

boundary  conditions  (LA.3),  the  variation  Jr  at  the  starting  and  finishing  points  5  =  0  and 
5  =  Z  +  JZ,  respectively,  should  be  zero,  i.e.  Jr(0)  =  0  and  Jr(Z  +  SI)  =  0.  The  last  condition 
is  used  to  calculate  the  variation  Jr  at  s  =  Z.  Namely,  from  Jr(Z  +  SI)  =  Jr(Z)  +  rSl  =  0  we 
obtain  Sr{l)  =  -f  JZ.  Using  the  last  equality,  the  variation  is  rearranged  in  the  form 

J  ■Sr  +  <l>5Xds 

4-  L  +  X<t>-  If  +  A  If  -r  ^^iSl 


Since  after  relaxing  constraint  (1.A.4),  all  three  variations  (Jx,  Sy,  Sz)  became  independent, 
the  necessary  conditions  for  an  extremum,  i.e.  S^  —  0,  are  reduced  to  the  constraint  (1.A.4) 
and  the  following  equations 


dL 

d  dL  d 

d 

11 

(1.A.6) 

dr 

ds  dr  ds 

dr 

L- 

dL,.d4> 
dr  dr 

•  r  =0. 

s=l 

(1.A.7) 

Vectorial  equation  (1.A.6)  has  the  first  integral.  Indeed,  the  scalar  product  of  (1.A.6)  with  f 
gives 

•r  =  0 

dr  ds  dr  ds  dr  ’ 

The  left-hand  side  of  this  equality  is  a  total  differential,  which  after  integration  becomes 


■  .  dL  .  .  d<l> 

L-r--^-X  r  -^ 
dr 


=  const. 


(1.A.8) 
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Lagrange  multiplier,  A(s),  is  derived  from  (1.A.8) 


AW- 


1/  -  r  •  +  Cl 

or 


.  d<t> 
dr 


(1A.9) 


where  cl  =  —const  is  an  unknown  constant.  Substitution  of  (1.A.9)  into  (1.A.6)  leads  to  the 
vectorial  differential  equation  for  determining  optimal  r 


dr 


±  ^+3- 

ds  dr  - 


— gr  L-t-^+cl 
dr 


=  0, 


(LA.10) 


which  along  with  the  constraint  (LA.4)  and  boundary  conditions  (1.A.3)  are  necessary  condi¬ 
tions  for  an  extremum.  Note  that  equations  (I.A.IO)  and  (1.A.4)  are  dependent  in  the  sense 
that  the  scalar  product  of  (I.A.IO)  with  f  is  reduced  to  r  •  =  0,  which  is  the  differential  of 

(1.A.4). 

In  the  case  when  constraint  I  <  U  is  active,  i.  e.  /  =  L,  equation  (1.A.7)  is  excluded  from 
determining  an  optimal  solution,  since  in  this  case  curve  total  length  is  fixed  and,  therefore,  the 
variation  SI  should  equal  zero  by  definition.  If  constraint  I  <  h  is  inactive,  then  from  (LA. 7) 
and  (1.A.8)  we  have  cl  =  0. 

However,  (I.A.IO)  with  (1.A.4)  and  (1.A.3)  are  only  the  necessary  conditions  for  an  optimal 
solution  to  solve  minimization  problem  (1.A.1)-(1.A.5),  since  (I.A.IO),  (1.A.4)  and  (1.A.3)  find 
an  extremal  trajectory,  which  either  minimizes  or  maximizes  the  functional. 

In  the  case  of  active  length  constraint  (variation  51  is  zero),  sufficient  condition  for  an  ex¬ 
tremal,  r*,  minimizing  functional  (1.A.2)  is  formulated 
for  all  r  sufficiently  close  to  r*  and 

(a)  all  r  sufficiently  close  to  r*  (weak  minimum); 

(b)  all  r  (strong  minimum); 
the  following  relation  holds 


u  u 

(L(r,  r)  -f  AV(f ))  ds  >  (L{r\  r*)  +  A^f"))  ds, 

0  0 

where  X*  is  given  by  (LA,9)  calculated  at  r*.  This  condition  is  reduced  to  verification  of  wether 
the  second  variation  of  the  functional  $(L,  A*,  I)  at  r*  is  greater  than  or  equal  to  zero. 
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1.  Introduction 

This  study  develops  a  general  approach  to  managing  risk  in  military 
applications  involving  stochasticity  and  uncertainties  in  distributions. 
Various  military  applications  such  as  intelligence,  surveillance,  planning, 
scheduling  etc.,  involve  decision  making  in  dynamic,  distributed,  and  im- 
certain  environments.  In  a  large  system,  multiple  sensors  may  provide 
incomplete,  conflicting,  or  overlapping  data.  Moreover,  some  compo¬ 
nents  or  sensors  may  degrade  or  become  completely  unavailable  (failures, 
weather  conditions,  battle  damage).  Uncertainties  in  combat  environ¬ 
ment  induce  different  kinds  of  risks  that  components,  sensors  or  armed 
units  are  exposed  to,  such  as  the  risk  to  be  damaged  or  destroyed,  risk 
of  mission  incompleteness  (e.g.,  missing  a  target)  or  failure,  risk  of  false 
target  attack  etc.  Therefore,  planning  and  operating  in  stochastic  and 
uncertain  conditions  of  a  modern  combat  require  robust  decision-making 
procedures.  Such  procedures  must  take  into  account  the  stochastic  na¬ 
ture  of  risk-inducing  factors,  and  generate  decisions  that  are  not  only  ef¬ 
fective  on  average  (in  other  words,  have  good  “expected”  performance), 
but  also  safe  enough  under  a  wide  range  of  possible  scenarios.  In  this 
regard,  risk  management  in  military  applications  is  similar  to  practices 
in  other  fields  such  as  finance,  nuclear  safety,  etc.,  where  decisions  tar¬ 
geted  only  at  achieving  the  maximal  expected  performance  may  lead  to 
an  excessive  risk  exposure.  However,  in  contrast  to  other  applications, 
distributions  of  the  stochastic  risk-inducing  factors  are  often  unknown 
or  uncertain  in  military  problems.  Uncertainty  in  distributions  of  risk 
parameters  may  be  caused  by  a  lack  of  data,  unreliability  of  data,  or 
the  specific  nature  of  a  risk  factor  (e.g.,  in  different  circumstances  a  risk 
factor  may  exhibit  different  stochastic  behavior).  Therefore,  decision 
making  in  military  applications  must  account  for  uncertainties  in  dis¬ 
tributions  of  stochastic  parameters  and  be  robust  with  respect  to  these 
uncertainties. 

In  this  project,  we  propose  a  general  methodology  for  managing  risk 
in  military  applications  involving  various  risk  factors  as  well  as  uncer¬ 
tainties  in  distributions.  We  build  our  approach  to  risk  management 
applications  on  the  CVaR  methodology,  which  is  a  relatively  new  devel¬ 
opment  (Rockafellar  and  Uryasev,  2000,  Rockafellar  and  Uryasev,  2001). 
The  approach  is  tested  with  several  stochastic  versions  of  the  Weapon- 
Target  Assignment  problem. 

The  report  is  organized  as  follows.  Section  2  develops  various  formula¬ 
tions  of  the  stochastic  Weapon-Target  Assignment  (WTA)  problem  with 
CVaR  constraints.  Results  of  numerical  experiments  for  one-stage  and 
two-stage  stochastic  WTA  problems  are  presented  in  Section  3.  The 
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Conclusions  section  summarizes  the  obtained  results  and  outlines  the 
directions  of  future  research.  The  Appendix  presents  key  theoretical  re¬ 
sults  on  risk  management  using  Conditional  Value-at-Risk  (CVaR)  risk 
measure,  and  describes  the  general  approach  to  controlling  risk  when 
distributions  of  risk  factors  are  uncertain. 

2.  Stochastic  Weapon- Target  Assignment 
Problem 

The  Weapon-Target  Assignment  (WTA)  problem  considers  the  opti¬ 
mal  assignment  of  weapons  to  targets  so  as  to  minimize  the  surviving 
value  of  targets.  The  WTA  problem  is  used  in  planning  environment 
that  features  a  whole  spectrum  of  uncertainties,  such  as  the  number  and 
types  of  targets  in  the  battle  space,  their  positions,  and  the  probability 
of  a  weapon  to  destroy  a  target  (e.g.  probability  of  kill).  To  gener¬ 
ate  robust  decisions,  one  must  account  for  these  uncertainties  and  the 
corresponding  risks.  In  this  section  we  present  two  formulations  of  the 
stochastic  Weapon-Target  Assignment  problem  that  address  the  uncer¬ 
tainties  in  a  weapon’s  probability  of  kill  and  in  the  number  of  targets. 

2.1.  Deterministic  WTA  Problem 

The  generic  formulation  of  the  Weapon  Target  Assignment  problem 
is  as  follows.  Given  the  set  of  targets  and  set  of  available  weapons,  one 
must  find  the  optimal  assignment  of  weapons  to  targets,  such  that,  for 
example,  the  damage  to  the  targets  is  maximized,  or  the  cost  of  the  op¬ 
eration  is  minimized.  The  WTA  formulation  that  maximizes  the  damage 
to  the  targets  (see,  for  example,  Manne,  1958,  denBroeger  et  al.,  1959, 
Murphey,  1999)  leads  to  a  non-linear  progranuning  problem  with  linear 
constraints  (NLP),  and  is  the  subject  of  future  research.  In  this  study 
we  adopt  another  setup,  where  the  total  cost  of  the  mission  (includ¬ 
ing  battle  damage  or  loss)  is  minimized,  while  satisfying  constraints  on 
mission  accomplishment  (i.e.,  destruction  of  all  targets  with  some  pre¬ 
scribed  probabilities).  We  assume  that  different  weapons  have  different 
costs  and  efficiencies,  and,  in  general,  each  may  have  a  “multishot”  ca¬ 
pacity  so  that  it  may  attack  more  than  one  target.  In  the  deterministic 
setup  of  the  problem  we  include  also  the  constraint  that  prescribes  how 
many  targets  a  single  weapon  can  attack. 

The  deterministic  WTA  problem  is 

^  I 

min  EE 

^  fc=l  t=l 
subject  to 
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'Y^Xik<mi,  i  =  l, ...,/, 

(lb) 

fc=l 

Xik<miVik,  i  =  l, ...,/,  k  —  l,...,K, 

(Ic) 

K 

(Id) 

k=l 

I 

1  -  J][(l >  difc,  =  (le) 

i=l 

^ik  ez+,  Vike{o,i}, 

where 

Xik  is  the  number  of  shots  to  be  fired  by  weapon  i  at  target  fc; 

Vik  =  1,  if  weapon  i  fires  at  target  k,  and  Vik  =  0  otherwise; 

Cjfc  is  the  cost  (including  the  battle  loss  or  damage)  of  firing  one  shot 
firom  weapon  i  at  target  fe;  Cfc  includes  the  relative  value  of  target 
k  with  respect  to  all  other  targets; 

rrii  is  the  shots  capacity  for  weapon  i\ 

ti  is  the  maximal  number  of  targets  which  can  be  attacked  by  weapon 

Pik  is  the  probability  of  destroying  target  k  by  firing  one  shot  firom 
weapon  i; 

dk  is  the  minimal  required  probability  for  destroying  target  k] 

Z  is  the  set  of  integer  niunbers,  and  Z"^  is  the  set  of  non-negative 
integers. 

The  objective  function  in  this  problem  equals  to  the  total  cost  of  the 
mission.  The  first  constraint,  (lb),  states  that  the  munitions  capacity 
of  weapon  i  cannot  be  exceeded.  The  second  and  the  third  constraints 
(Ic)  and  (Id)  are  responsible  for  not  allowing  weapon  i  to  attack  more 
than  ti  targets,  where  U  <  K.  The  last  constraint  (le)  ensures  that  after 
all  weapons  are  assigned,  target  k  is  destroyed  with  probability  not  less 
than  dk-  Note  that  this  non-linear  constraint  can  be  linearized: 

I 

-  Pik)  Xik  -  in(l  -  dk)  <  0.  (2) 

t=i 

In  this  way  the  deterministic  WTA  problem  (la)  can  be  formulated  as 
a  linear  integer  programming  (IP)  problem. 
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2.2.  One-Stage  Stochastic  WTA  Problem  with 
CVaR  Constraints 

In  reaHife  situations  many  of  the  parameters  in  model  (la)-(le)  are 
not  deterministic,  but  stochastic  values.  For  example,  the  probabilities 
Pik  of  destroying  target  k  may  depend  upon  battle  situation,  weather 
conditions,  and  so  on,  and  consequently,  may  be  treated  as  being  uncer¬ 
tain.  Similarly,  the  cost  of  firing  c^fc,  which  includes  battle  loss/damage, 
may  also  be  a  stochastic  parameter.  The  number  of  targets  K  may  be 
uncertain  as  well. 

First,  we  consider  a  one-stage  Stochastic  Weapon-Target  Assignment 
(SWTA)  problem,  where  the  uncertainty  is  introduced  into  the  model 
by  assuming  that  probabilities  pik  are  stochastic  and  dependent  on  some 
random  parameter 


Pik  —  Pik{^)- 


In  accordance  to  the  described  methodology  of  managing  uncertainties 
and  risks  in  military  applications,  we  model  the  stochastic  behavior  of 
probabilities  pifc  using  scenarios.  Namely,  probabilities  Pifc(0  take  differ¬ 
ent  values  Piki^s)  =  Piks^  «  =  -j  S  under  S  different  scenarios.  Such  a 
scenario  set  may  be  constructed,  for  example,  by  utilizing  the  historical 
observations  of  weapons’  efficiency  in  different  environments,  or  by  using 
simulated  data,  experts’  opinions  etc. 

To  control  risks  we  use  Conditional  Value-at-Risk  (CVaR)  approach. 
A  general  risk  management  approach  with  CVaR  functions  is  described 
in  Appendix,  We  replaced  the  last  constraint  in  (la)  by  a  CVaR  con¬ 
straint,  where  the  loss  function  takes  a  positive  value  if  the  probability 
of  destroying  target  k  is  less  than  dk' 

Lk{x,  0  =  ^ik  -  ln(l  -  4),  (3) 

and  takes  a  negative  value  otherwise.  The  CVaR  constraint  with  confi¬ 
dence  level  a  bounds  the  (weighted)  average  of  (1  —  a)  •  100%  highest 
losses.  In  our  case,  allowing  small  positive  values  of  loss  function  (3)  for 
some  scenarios  implies  that  for  these  scenarios  target  k  is  destroyed  with 
probability  slightly  less  than  4?  which  may  still  be  acceptable  from  a 
practical  point  of  view. 

Except  for  the  constraint  on  the  target  destruction  probability,  the 
one-stage  Stochastic  WTA  problem  is  identical  to  its  deterministic  pre- 
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decessor: 

K  I 

min  EE  (4) 

^  k=l  2=1 
subject  to 
K 

'Y^Xik<mu  i  =  l, 

fc=i 

^ih  —  i  —  1)  •••)  ^  1)  •••) 

fc=l 

CVaRa  [Lfe(x,  0]  <  C'fc.  =  1.  ^• 

Here  a  is  the  confidence  level,  Ck  are  some  (small)  constants,  and  all 
other  variables  and  parameters  are  defined  as  before.  As  demonstrated 
in  (11),  for  the  adopted  scenario  model  with  probabilities  pik,  the  CVaR 
constraint  for  the  fc-th  target 

CVaRa[Ljk(a;,e)]<C'fc 

is  represented  by  a  set  of  linear  inequalities: 

I 

- piks) Xik  -  ln{l  -  dk)  -  Ck  < ‘^sk,  s  =  l,...,S, 

2=1 

S 

a + (1  -  ockr^s-'^Y^wsk  <  Ck,  (5) 

5=1 

Ck^T^,  Wsk>0,  s  =  l,...,S,  k  =  l,...,K. 

Thus,  the  one-stage  Stochastic  WTA  problem  can  be  formulated  as  a 
mixed-integer  programming  (MIP)  problem: 

K  I 

min  EE  ^k  ^ik  (6) 

^  fe=l  2=1 
subject  to 
K 

^  ^  ^ik  ^  ^2}  ^  ^  •••1 

A;=l 

^ik  ^  ^2  ^2^:5  i  =  Ij  /, 

K 

^  ^  '^ik  ^  ^2)  ^  ^  1? 

k=l 


k  = 
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I 

^  lii(l  -  Piks)  Xik  -  ln(l  -  dk)  -Ck<  ‘Wsk, 
s  =  1, k  =  1,...,^^", 

Ca:  +  (1  “  ^  '^sk  ^  ?  A:  =  l,...,i^5 

s=l 

Xifc  e  2:+,  Vik  €  {0, 1},  a  ^  >  0, 

5  =  1,. ..,5,  i  =  i, = 

Note  that  different  values  of  probability  Pik  represent  the  uncertainty 
in  the  distributions  of  stochastic  parameters  discussed  in  the  previous 
section.  Indeed,  different  values  of  probability  pik  imply  different  prob¬ 
ability  measures  for  the  random  variable  associated  with  the  event  of 
destroying  target  k  by  firing  one  unit  of  munitions  by  weapon  i.  In  ef¬ 
fect,  CVaR  constraint  (4)  is  a  risk  constraint  that  incorporates  multiple 
probability  measures. 

2.3.  Two-Stage  Stochastic  WTA  Problem  with 
CVaR  constraints 

In  this  section  we  consider  a  more  complex,  but  also  more  realistic 
two-stage  Stochastic  WTA  problem,  where  the  uncertain  parameter  is 
the  number  of  targets  to  be  destroyed. 

This  problem  is  more  realistic  since  it  models  the  effect  of  target 
discovery  as  being  dynamic;  that  is,  not  all  targets  are  known  at  any 
single  instance  of  time.  To  address  this  type  of  uncertainty,  we  need  to 
modify  our  notation. 

Consider  I  weapons  are  deployed  in  some  bounded  region  of  interest 
and  interval  of  time  T  with  the  goal  of  finding  targets  and  then,  once 
found,  attacking  those  targets.  If  we  delay  all  assignments  of  weapon 
shots  until  to  targets  until  the  final  time  T,  then  we  have  a  deterministic, 
“static”  WTA  problem  as  in  (la)-(le).  If,  on  the  other  hand,  we  assume 
that  weapons  have  at  least  2  opportunities  to  shoot  during  the  interval 
r,  then  the  WTA  problem  is  dynamic.  In  the  later  case  we  have  the 
opportunity  to  avoid  expending  all  our  shots  at  targets  discovered  early 
in  T  by  explicitly  modeling  the  number  of  imdiscovered  targets  in  the 
objective  function. 

Assume  that  K  now  represents  the  number  of  categories  of  targets 
(the  targets  may  be  categorized,  for  example,  by  their  importance,  vul¬ 
nerability,  etc). 

We  will  assume  the  problem  has  2  stages.  That  is,  at  any  given  point  in 
time,  we  may  always  partition  all  targets  into  those  thus  far  determined 
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and  those  that  we  conjecture  to  exist  but  have  not  yet  found.  Our 
conjecture  may  be  based  on  evidence  obtained  by  prior  reconnaissance 
of  the  region  of  interest.  At  some  arbitrary  time  0  <  r  <  T  assume  that 
there  are  rifc  detected  targets  and  r}k  undetected  targets  in  each  category 
k  =  1) K.  Thus  we  have  two  clearly  identified  stages  in  our  problem: 
in  the  first  stage  one  has  to  destroy  the  targets  known  at  time  r ,  in  the 
second  stage  one  must  destroy  the  targets  that  we  conjecture  will  be 
found  by  time  T.  In  other  words,  one  needs  to  make  an  assignment  of 
weapons  that  will  allow  for  the  destruction  of  the  targets  known  at  time 
T  while  reserving  enough  munition  capacity  for  destroying  the  targets 
we  expect  to  find  in  r  <  t  <  T. 

Setup  of  the  two-stage  stochastic  WTA  problem  can  be  considered  as  a 
part  of  a  moving  horizon  or  quasi-multistage  stochastic  WTA  algorithm, 
where  the  WTA  problem  with  many  time  periods  is  solved  by  recursive 
application  of  a  two-stage  algorithm  (Murphey,  1999). 

To  simplify  the  problem  setup,  we  remove  the  constraint  on  the  num¬ 
ber  of  targets  a  single  weapon  can  attack  (the  second  and  third  con¬ 
straints  in  problems  (la)),  since  this  constraint  makes  the  problem  much 
too  combinatorial.  Also,  we  assume  that  the  probabilities  pik  are  known 
{not  random),  so  that  the  only  stochastic  parameters  in  the  two-stage 
SWTA  problem  are  the  numbers  of  undetected  (second-stage)  targets 
77fc)  A:  =  1,  ...j  K. 

We  model  the  uncertainty  in  the  number  of  targets  at  the  second  stage, 
by  we  introducing  a  scenario  model,  where  under  scenario  s  €  {1, ...,  5} 
there  are  77^(5)  =  rjks  undetected  targets  in  category  k. 

The  first-  and  second-stage  decision  variables  are  defined  as  follows: 

Xik  is  the  number  of  munitions  to  be  fired  by  weapon  i  at  a  single 
target  in  category  k  during  the  first  stage; 

yik{s)  is  the  number  of  munitions  to  be  fired  by  weapon  i  at  a  single 
target  in  category  k  during  the  second  stage  scenario  s. 

Note  that  the  same  decision  is  made  for  all  targets  within  a  category, 
i.e.,  once  weapon  i  fires,  say,  2  missiles  at  a  specific  target  in  category 
fe,  it  must  fire  2  missiles  at  every  other  target  in  this  category. 

The  recursive  formulation  of  the  two-stage  stochastic  WTA  problem 
is 

min  |y^ nkCjkXjk  +  Erj[Q{x, 77)]  1  (7a) 

U=li=l  ) 

subject  to 
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K 

fe=l 

/ 

ln(l  -  Pik)  Xik  -  ln(l  -  dk)  <  Sik,  A:  =  1, K,  (7c) 

i=l 

j2^ik<C,  (7d) 

fc=i 

^Ik  —  ^  1)  •••) 

where  the  recourse  function  Q{x,  rf)  is  the  solution  of  the  problem 

{  K  I  I  \ 

Q{x,  Tj)  =  min 

*'  U=li=l  i=l  ) 

subject  to 
K 

{rikXik  +  mi^)  yik{s))  <mi  +  Si,  V i,  (8b) 

A;=l 

I 

ln(l  -  Pik)  yik{s)  -  ln(l  -  4)  -Ck<  '^k{s),  V  fc,  s,  (8c) 

i=l 

S 

Ck  +  (}-0Ck)^^S^^'^Wk{s)  <e2k,  vfe,  (8d) 

s=l 

K 

^(£ifc  +  £2fc)  ^  C*,  (8e) 

fc=i 

yifc(s),  €  -2+,  Wk{s),  e2k  >0,  Cfc  €  7?-,  M  »  1. 


Let  us  discuss  the  recourse  problem  (7a)-(8e).  As  before,  we  minimize 
the  total  cost  of  the  mission.  The  first  constraint  (7b)  is  the  munitions 
capacity  constraint.  The  second  constraint,  (7c),  allows  a  first-stage 
target  in  category  k  to  siurvive  with  (small)  error  eifc,  and  the  third 
constraint  (7d)  bounds  the  sum  of  errors  Sifc  by  some  (small)  constant 
C. 

In  the  recourse  function  (8a)  the  first  constraint  (8b)  requires  the 
weapon  i  to  not  exceed  its  munitions  capacity  while  destroying  the  first- 
and  second-stage  targets.  The  possible  infeasibility  of  the  munitions 
capacity  constraint  can  be  relaxed  using  auxiliary  variables  5i  that  enter 
the  objective  function  with  cost  coefficient  M  The  second  and  third 
constraints  (8c)-(8d)  form  a  CVaR  constraint  that  controls  the  failure 
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of  destroying  second-stage  targets  with  the  prescribed  probabilities  dfc. 
Similarly  to  the  deterministic  constraint  in  (7a),  CVaR  of  failure  to 
destroy  a  second-stage  target  in  category  k  is  bounded  by  (small)  error 
variable  e2k-  The  total  sum  of  errors  eifc  and  e2fc  at  both  stages  is 
bounded  by  small  constant  C,  which  makes  possible  a  tradeoff  between 
the  degree  of  mission  accomplishment  at  the  first  and  second  stages. 

The  extensive  form  of  the  two-stage  SWTA  problem  (7a)— (8a)  is 


{K  I  ^  S  K  J  I 

^'^nkCikXik  + 

fc=l  t=l  ^  s=l  fe=l  i=l  i=l 

subject  to 
K 

(nfc  Xik  +  Vks  Viki^))  ^  V s, 

*=1 

I 

E‘“(i  -  Pik)  Xik  -  -  dk)  <  £ik, 

i=l 
I 

ln(l  -  Pik)  Vikis)  -  ln(l  -  dk)  -  Cifc  < 

i=l 

S 

Cfc  -b  (1  -  ak)~^S~^  Wks  <  £2k,  V  k, 

S=1 

K 

y~i(gifc + £2fc)  ^  C) 

k=l 

Xik,  yik{s),  Si  €  '2'^)  ^fes.  eU)  ^2k  Cfc  €  T?.,  M  »  1. 


(9) 


The  two-stage  stochastic  WTA  problem  is  also  a  MIP  problem. 

3.  Numerical  results 

In  this  section  we  present  and  discuss  numerical  results  obtained  for 
both  one-stage  and  two-stage  stochastic  WTA  problems.  The  algorithms 
for  solving  deterministic,  one-  and  two-stage  stochastic  WTA  problems 
were  implemented  in  C-t-+,  and  we  used  CPLEX  7.0  Callable  Library  to 
solve  the  corresponding  IP  and  MIP  problems.  We  used  simulated  data 
(sets  of  weapons  and  targets,  the  corresponding  costs  and  probabilities 
etc.)  for  testing  the  implemented  algorithms. 
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3.1.  Single-stage  deterministic  and  stochastic 
WTA  problems 

For  the  deterministic  and  one-stage  stochastic  WTA  problems  we  used 
the  following  data: 

■  5  targets  {K  =  5) 

■  5  weapons,  each  with  4  shots  (/  =  5,  mj  =  4) 

■  any  weapon  can  attack  any  target  {U  =  5), 

■  probabilities  put  and  costs  Cik  depend  only  on  the  weapon  index  i: 

Pik  ~  Pi  5  ^ 

■  all  targets  have  to  be  destroyed  with  at  least  probability  95%  {dk  = 
0.95) 

■  the  confidence  levels  in  CVaR  constraint  are  0.90 

■  there  are  20  scenarios  (S  =  20)  for  probabilities  Pik{s)  in  the  one- 
stage  SWTA  problem;  all  scenarios  are  equally  probable. 

According  to  the  aforementioned,  we  used  simulated  data  for  probabil¬ 
ities  piks  and  costs  Cik^  It  was  assumed  that  probabilities  piks  =  Pis  are 
imiformly  distributed  random  variables,  and  the  Fig.  3.1  displays  the 
relation  between  the  cost  of  missile  of  weapon  i  and  its  efficiency  (i.e., 
probability  to  destroy  a  target): 


Figure  1.1.  Dependence  between  the  cost  and  efficiency  for  different  types  of  weapons 
in  one-stage  SWTA  problem  (6)  deterministic  WTA  problem  (la). 


On  this  graph,  diamonds  represent  the  average  probability  of  destroying 
a  target  by  firing  one  shot  fi:om  weapon  t,  and  the  horizontal  segments 
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represent  the  support  for  random  variable  Pik{0  =  The  average 

probabilities 

1  ^ 

Pik  —  "g  y^.Piks 

5=1 


were  used  for  pik  in  the  deterministic  problem  (la). 

The  efficiency  and  cost  of  weapons  1  to  5  increase  with  the  index 
of  weapon,  i.e.,  Weapon  1  is  the  least  efficient  and  cheapest,  whereas 
Weapon  5  is  the  most  precise,  but  also  most  expensive  one. 

Tables  1.1  and  1.2  represent  the  optimal  solutions  (variables  Xik)  of 
the  deterministic  and  one-stage  stochastic  WTA  problems. 


Table  1.1.  Optimal  solution  of  the  deterministic  WTA  problem  (la) 


Target 

T1 

T2 

T3 

T4 

T5 

Total  shots 

Weapon  1 

0 

2 

1 

0 

1 

4 

Weapon  2 

0 

1 

2 

0 

0 

3 

Weapon  3 

1 

0 

0 

1 

1 

3 

Weapon  4 

1 

0 

0 

1 

1 

3 

Weapon  5 

0 

0 

0 

0 

0 

0 

Table  1.2.  Optimal  solution  of  the  one-stage  stochastic  WTA  problem  (4),  (6) 

Target  T1 

T2 

T3 

T4 

T5 

Total  shots 

Weapon  1  0 

1 

1 

0 

1 

3 

Weapon  2  0 

0 

1 

1 

1 

3 

Weapon  3  2 

0 

0 

1 

0 

3 

Weapon  4  0 

1 

1 

0 

1 

3 

Weapon  5  1 

1 

0 

1 

0 

3 

One  can  observe  the  difference  in  the  solutions  produced  by  determin¬ 
istic  and  stochastic  WTA  problems:  the  deterministic  solution  does  not 
use  the  most  expensive  and  most  precise  Weapon  5,  whereas  the  stochas¬ 
tic  solution  of  problem  (6)  with  CVaR  constraint  uses  this  weapon.  It 
means  that  the  CVaR-constrained  solution  of  problem  (6)  represents  a 
more  expensive  but  safer  decision. 

On  a  different  dataset,  we  obtained  a  similar  result:  the  optimal  so¬ 
lution  of  the  stochastic  problem  with  CVaR  constraints  did  not  use  the 
cheapest  and  the  most  unreliable  weapon,  whereas  the  deterministic  so¬ 
lution  used  it. 
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We  have  also  performed  testing  of  the  deterministic  solution  under 
different  scenarios.  The  deterministic  solution  failed  to  destroy  more 
than  one  target  under  13  of  20  scenarios. 

This  example  highlights  the  importance  of  using  risk  management  pro¬ 
cedures  in  military  decision-making  applications  involving  uncertainties. 

3.2.  Two-Stage  Stochastic  WTA  Problem 

For  the  two-stage  stochastic  WTA  problems  we  used  the  following 
data: 

■  3  categories  of  targets  {K  =  3) 

■  4  weapons,  each  with  15  shots  (/  =  4,  mj  =  15) 

■  probabilities  pik  and  costs  Cik  depend  only  on  the  weapon  index  i: 
Pik  ~  Pi'i  ^  ^ 

■  all  targets  have  to  be  destroyed  with  probability  95%  (dfc  =  0.95) 

■  the  confidence  levels  in  CVaR  constraint  are  equal  0.90 

■  there  are  15  scenarios  {S  =  15)  for  the  number  of  imdetected  tar¬ 
gets  r]ks  (for  each  k,  the  number  of  undetected  targets  rjks  is  a 
random  integer  between  0  and  5);  all  scenarios  are  equally  proba¬ 
ble. 


Figure  1.2.  Dependence  between  the  cost  and  efficiency  for  different  types  of  weapons 
in  two-stage  SWTA  problem  (9). 


For  the  probabilities  pik  in  the  two-stage  problem,  we  used  the  first 
four  average  probabilities  from  the  deterministic  WTA  problem,  and 
the  efficiency-cost  dependence  is  shown  in  Fig.  3.2. 
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Tables  1.3  to  1.5  illustrate  the  optimal  solution  of  the  problem  (9). 
Table  1.3  contains  the  first-stage  decision  variables  Xik^  and  Tables  1.4 
and  1.5  display  the  second-stage  variables  yikis)  for  scenarios  s  =  1  and 
s  =  2,  just  for  illustrative  purposes. 

Similarly  to  the  analysis  of  the  one-stage  stochastic  WTA  problem,  we 
compared  the  scenario-based  solution  of  problem  (9)  with  the  solution 
of  the  “deterministic  two-stage”  problem,  where  the  number  of  second- 
stage  targets  in  each  category  is  taken  as  the  average  over  15  scenarios. 
The  comparison  shows  that  the  solution  based  on  the  expected  infor¬ 
mation  leads  to  significant  munitions  shortages  in  5  of  15  (i.e.,  33%) 
scenarios,  and  consequently  to  failing  the  mission  at  the  second  stage. 
Recall  from  the  analysis  of  the  one-stage  SWTA  problem  that  the  solu¬ 
tion  based  on  the  expected  information  also  exhibited  poor  robustness 
with  respect  to  different  scenarios.  Indeed,  solutions  that  use  only  the 
expected  information,  are  supposed  to  perform  well  on  average,  or  in 
the  long  run.  However,  in  military  applications  there  is  no  long  run,  and 
therefore  such  solutions  may  not  be  robust  with  respect  to  many  possible 
scenarios. 


Table  1,3.  First-stage  optimal  solution  of  the  two-stage  stochastic  WTA  problem 


Category 

K1 

K2 

K3 

#  of  detected  targets 

3 

5 

2 

Weapon  1 

0 

0 

0 

Weapon  2 

0 

0 

0 

Weapon  3 

1 

1 

1 

Weapon  4 

1 

1 

1 

Table  1.4-  First-stage  optimal  solution  of  the  two-stage  stochastic  WTA  problem  (9) 
for  the  first  scenario 


Category 

K1 

K2 

K3 

#  of  undetected  targets 

1 

4 

2 

Weapon  1 

0 

0 

2 

Weapon  2 

0 

0 

1 

Weapon  3 

1 

1 

0 

Weapon  4 

1 

1 

0 
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Table  1.5.  Second-stage  optimal  solution  of  the  two-stage  stochastic  WTA  problem 
(9)  for  the  second  scenario 


Category 

K1 

K2 

K3 

#  undetected  of  targets 

3 

5 

3 

Weapon  1 

2 

0 

2 

Weapon  2 

1 

0 

1 

Weapon  3 

0 

1 

0 

Weapon  4 

0 

1 

0 

Thus,  solutions  of  both  one-stage  and  two-stage  SWTA  problems  con¬ 
firm  the  general  conjecture  on  the  potential  importance  of  exploiting 
stochastic  models  and  risk  management  in  military  applications. 

4.  Conclusions 

We  have  presented  an  approach  to  managing  risk  in  stochastic  en¬ 
vironments,  where  distributions  of  stochastic  parameters  are  imcertain. 
This  approach  is  based  on  the  methodology  of  risk  management  with 
Conditional  Value-at-Risk  risk  measure  developed  by  Rockafellar  and 
Uryasev,  2000,  2001.  Although  the  presented  approach  has  been  used 
to  solve  one-stage  and  two-stage  stochastic  Weapon-Target  Assignment 
problems,  it  is  quite  general  and  can  be  applied  to  wide  class  of  problems 
with  risks  and  uncertainties  in  distributions.  Among  the  directions  of 
future  research  we  emphasize  consideration  of  a  stochastic  WTA  prob¬ 
lem  in  NLP  formulation,  where  the  damage  to  the  targets  is  maximized 
while  constraining  the  risk  of  false  target  attack. 
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5.  Appendix.  The  General  Approach  to  Risk 
Management 

Presence  of  uncertainty  in  a  decision-making  model  leads  to  the  prob¬ 
lem  of  estimation  and  managing/controlling  of  risk  associated  with  the 
stochastic  parameters  in  the  model.  Over  the  recent  years,  risk  manage¬ 
ment  has  evolved  into  a  sophisticated  discipline  combining  both  rigorous 
and  elegant  theoretical  results  and  practical  effectiveness  (this  especially 
applies  to  the  risk  management  in  finance  industry).  Generally  speaking, 
risk  management  is  a  set  of  activities  aimed  at  reducing  or  preventing 
high  losses  incurred  firom  an  incorrect  decision.  The  losses  (e.g.,  dam¬ 
ages,  failures)  in  a  system  are  quantified  by  a  loss  function  L{x,  that 
depends  upon  decision  vector  x  and  a  stochastic  vector  ^  standing  for 
uncertainties  in  the  model.  Assuming  for  now  that  a  distribution  of  the 
parameter  ^  is  known,  it  is  possible  to  determine  the  distribution  of  the 
loss  function  L{x,C)  (see  Fig.  1.3).  Then,  the  problem  of  preventing 
high  losses  is  a  problem  of  controlling  and  shaping  the  loss  distribution 
and,  more  specifically,  its  right  tail,  where  the  high  losses  reside.  To 
estimate  and  quantify  the  losses  in  the  tail  of  the  loss  distribution,  a  risk 
measure  has  to  be  specified.  In  particular,  a  risk  measure  introduces  the 
ordering  relationships  for  risks,  so  that  one  is  able  to  discriminate  “less 
risky”  decisions  from  the  “more  risky”  ones^.  The  appropriate  choice  of 
a  risk  measure  is,  in  most  cases,  dictated  by  the  nature  of  imcertainties 
and  risks  in  the  problem  at  hand.  In  military  applications,  for  example, 
one  usually  deals  with  the  probabilities  of  events,  such  as  the  proba¬ 
bility  to  hit  a  target,  the  probability  to  detect  the  enemy’s  aircraft. 


‘Artzner  et  al.,  1999,  have  introduced  a  concept  of  “ideal”,  or  coherent,  risk  measure.  A 
coherent  risk  measure,  which  satisfies  to  a  set  of  axioms  developed  in  this  paper,  is  expected 
to  produce  “proper*’  and  “consistent”  estimates  of  risk. 


REFERENCES 


79 


and  so  on.  Therefore  percentile  risk  measures  that  represent  the  risk 
in  terms  of  percentiles  of  the  loss  distribution  are  particularly  suitable 
for  the  risk  management  in  military  applications.  Popular  percentile  risk 
measures  include  Value- at-Risk  (VaR),  Conditional  Value-at-Risk,  Max¬ 
imum  Loss,  and  Expected  Shortfall.  Figure  1.3  displays  some  of  these 
measures;  Value-at-Risk  with  confidence  level  a  (a-VaR),  which  is  the 
a-percentile  of  loss  distribution,  Maximum  Loss  (‘1.0-percentile”  of  loss 
distribution),  and  a-CVaR,  which  may  be  thought  of  as  the  expectation 
of  losses  exceeding  a- VaR. 


Figure  1.3.  Loss  function  distribution  and  different  risk  measures. 


We  build  our  approach  to  risk  management  in  military  applications 
on  the  CVaR  methodology,  which  is  a  relatively  new  development  ( 
Rockafellar  and  Uryasev,  2000,  Rockafellar  and  Uryasev,  2001).  This 
section  presents  the  general  framework  of  risk  management  using  Con¬ 
ditional  Value-at-Risk,  and  extends  it  to  the  case  when  the  distributions 
of  stochastic  parameters  are  not  certain. 

5.1.  Risk  Management  Using  Conditional 
Value-at-Risk 

Consider  a  loss  function  L(®,  depending  on  a  decision  vector  x  and 
a  stochastic  vector  and  its  cumulative  distribution  function  (c.d.f.) 


nx,0  =  P[L{x,0<Cl 
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Then  the  a- VaR  (Value-at-Risk  at  confidence  level  a)  function  Ca(®) 
corresponding  to  loss  L(x,  is 


Ca(a:)  =  C)  >  «}• 

An  “intuitive”  definition  of  Conditional  Value-at-Risk  with  confidence 
level  a  (o-CVaR)  presents  it  as  conditional  expectation  of  losses  ex¬ 
ceeding  the  a- VaR  level.  This  definition  is  correct,  however,  only  for 
continuously  distributed  loss  functions.  For  loss  functions  with  gen¬ 
eral  non-continuous  distributions  the  a-CVaR  function  ^a{x)  is  defined 
as  the  expected  value  of  random  variable  Za  (Rockafellar  and  Uryasev, 
2001): 

Mx)  =  CVaRa[T(a:,^)]  =  EM, 

where  c.d.f.  C)  of  Za  has  the  form 


[i(x,C)  -  a]/[l  -n], 


C  < 

C  >  Ca(®)- 


In  (Rockafellar  and  Uryasev,  2001),  it  was  shown  that  n-CVaR  can  be 
expressed  as  a  convex  combination  of  a- VaR  and  conditional  expectation 
of  losses  strictly  exceeding  a- VaR: 


<f>a{x)  =  Aa(x)  Ca(a;)  +  [1  -  Aa(x)]  <?!>J(x),  (10) 


where 

^+(x)  =  E[Z(x,0  I  L{x,0  >  Ca(x)],  (11) 

which  is  also  known  as  “upper  CVaR”  or  Expected  Shortfall,  and 

Aa(x)  =  [’5'(x,  Ca(®))  -  0:]/[l  -  0  ^  -^0(3:)  <  1. 

Similar  to  (11),  another  percentile  risk  measure,  called  “lower  CVaR”, 
or  CVaR“,  can  be  defined: 

<t>-{x)  =  E[L(x,0  I  Lix,0  >  Ux)]- 

Then,  as  it  was  shown  in  (Rockafellar  and  Uryasev,  2001),  the  introduced 
risk  functions  satisfy  the  following  inequality: 

Ca(x)  <  <p~{x)  <  4>a{x)  <  (l>iix). 

The  conditional  Value-at-Risk  function  <^a{x)  has  the  following  prop¬ 
erties  (Rockafellar  and  Uryasev,  2000,  Rockafellar  and  Uryasev,  2001, 
Acerbi  and  Tasche,  2001): 
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■  CVaR  is  continuous  with  respect  to  confidence  level  a  (VaR, 
CVaR''",  CVaR“  may  be  discontinuous  in  a); 

■  CVaR  is  convex  in  a  and  x,  provided  that  the  loss  function  L{x, 

is  convex  in  x  (VaR,  CVaR+,  CVaR"  are  generally  non-convex  in 

x); 

■  CVaR  is  coherent  in  the  sense  of  Artzner  et  al.,  1999; 

■  unlike  VaR,  CVaR  has  stable  statistical  estimates; 

■  in  the  case  of  a  continuous  loss  distribution  CVaR  coincides  with 
CVaR'*'  and  CVaR",  and  represents  the  conditional  expectation  of 
losses  exceeding  VaR.  ; 

From  the  viewpoint  of  managing  and  controlling  of  risk,  the  most  impor¬ 
tant  property  of  CVaR,  which  distinguishes  it  from  all  other  percentile 
risk  measures,  is  the  convexity  with  respect  to  decision  variables,  which 
permits  the  use  of  convex  programming  for  minimizing  CVaR.  If  the  loss 
function  L(x,  can  be  approximated  by  a  piecewise  linear  function,  the 
procedure  of  controlling  or  optimization  of  CVaR  is  reduced  to  solving 
Linear  Programming  (LP)  problem. 

The  techniques  for  optimizing  CVaR  when  the  loss  distribution  is 
discrete  are  of  special  importance  for  mihtary  apphcations,  as  will  be 
demonstrated  in  Section  3. 

Assume  that  there  are  S  possible  realizations  (scenarios)  ...,^s  of 
vector  ^  with  probabilities  tt*  (Z)f=i  tTs  =  1),  then  in  the  optimization 
problem  with  multiple  CVaR  constraints 

max  g(x) 
subject  to 

CVaRajL(rr,0]<C'n,  n  = 

where  g{x)  is  some  performance  function  and  X  is  a  convex  set,  each 
CVaR  constraint  may  be  replaced  by  a  set  of  inequalities 

Cn  ^  S  =  1)  •••)  Sj 

S 

Cn  +  (1  ^n)  ^  '^s'^ns  ^  (12) 

s=l 

(^fi  €  7?^,  ^  72."^,  S  =  1, S, 

where  TZ  is  the  set  of  real  numbers,  and  is  the  set  of  non-negative 
real  numbers,  and  Wns  sire  auxiliary  variables.  If  in  the  optimal  solution 
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the  n-th  CVaR  constraint  is  active,  then  the  corresponding  variable  Cn 
is  equal  to  an- VaR  (i.e.,  an- th  percentile  of  the  loss  distribution). 

We  also  note  that  in  the  case  when  the  behavior  of  stochastic  param¬ 
eter  ^  can  be  represented  by  a  scenario  model  s  =  1, 5}  with 
equally  probable  scenarios  (tts  =  the  concept  of  CVaR  acquires 

especially  simple  and  transparent  interpretation.  Namely,  if  (for  a  fixed 
x)  the  scenarios  ...,  ^5  are  indexed  such  that  L(x,  ^1)  <  ■■■  <  L{x,  ^s), 
then  a-CVaR  equals  the  weighted  average  of  losses  for  the  [(1  —  a) S'] 
worst  scenarios: 


(1  _  a)  5 

where  number  Sa  is  such  that 

5  —  5a  <  (1  Ol)S  <  S  —  Sa  +  l- 

In  the  risk  management  methodology  discussed  above  the  distribution 
of  stochastic  parameter  ^  is  considered  to  be  known.  The  next  subsection 
extends  the  presented  approach  to  the  case,  when  the  distribution  of 
stochastic  parameters  in  the  model  is  not  certain. 

5.2.  Risk  Management  Using  CVaR  in  the 

Presence  of  Uncertainties  in  Distributions 

The  general  approach  to  managing  risks  in  an  uncertain  environment, 
where  the  distributions  of  stochastic  parameters  are  not  known  for  sure, 
can  be  described  as  follows.  Suppose  that  we  have  some  performance 
function  F{x,^),  dependent  on  the  decision  vector  x  E  X  and  some 
random  vector  ^  €  S,  whose  distribution  is  not  known  for  certain.  We 
assume  that  the  actual  realization  of  vector  ^  may  come  from  different 
distributions  ©1,..., 0jv-  The  vector  ^  stands  for  the  uncertainties  in 
data  that  make  it  impossible  to  evaluate  the  efficiency  F{x,^)  of  the 
decision  for  sure.  Thus,  there  always  exists  a  possibility  of  making  an 
incorrect  decision,  and,  consequently,  suffering  loss,  damage,  or  failing 
the  mission.  If  the  loss  in  the  system  is  evaluated  by  function  L(x,^), 
then  risk  of  high  losses  can  be  controlled  using  CVaR  constrains.  Let 
formulate  the  problem  of  maximizing  the  expected  performance  function 
F{x,  subject  to  some  operational  constraints  Ax  <  b  and  CVaR  risk 
constraints.  Due  to  the  unknown  distribution  of  vector  we  are  unable 
to  find  the  expectation  Ee[F(x,  ^)].  Therefore,  being  on  the  conservative 
side,  we  want  the  decision  x  to  be  optimal  with  respect  to  each  measure 


{Sa  ~  CxS)  L(x,^Sa)  ^  > 

S=Sa+l 
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©n,  and  this  leads  to  the  following  maxi-min  problem: 

x€X  ©n, 

subject  to 
Ax  <  b, 

CyaRa[L{x,0  I  ©n]  <  C,  n  =  1, N, 


(13) 


where  multiple  CVaR  constraints  with  respect  to  different  measures 
control  the  risk  for  high  losses  L{x^^)  to  exceed  some  threshold  C.  In 
formulation  (13)  we  assume  that  the  performance  function  F  is  concave 
in  and  the  loss  function  L  is  convex  in  x.  These  assumptions  are 
not  restrictive;  on  the  contrary,  they  indicate  that  given  more  than  one 
decision  with  equal  performance  one  favors  safer  decisions  over  the  riskier 
ones. 

Model  (13)  explains  how  to  handle  the  risk  of  generating  an  incorrect 
decision  in  an  uncertain  environment.  In  military  applications,  different 
types  of  risks  and  losses  may  be  explicitly  involved,  for  example,  along 
with  loss  function  L(a;,^)  one  may  consider  a  loss  function  R{x^  ^)  for 
the  risk  of  false  target  attack.  Control  for  this  type  of  risk  can  also  be 
included  in  the  model  by  a  similar  set  of  CVaR  constraints: 


max  min  E©^[F(x,0] 

2c€X  ©n, 

subject  to 

CVaRai  [X(x,  0  I  e„]  <  Cl,  n  =  1, N, 

CVaRaARi^,  0 1  ©n]  <  C2,  n  =  1, N. 

In  the  next  sections  we  test  the  presented  approach  to  risk  management 
in  military  applications  on  the  Weapon-Target  Assignment  problem. 
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1.  Introduction 

Let  G{V,  E)  be  an  undirected  graph  with  set  of  vertexes  F  =  {1, 2, . . . ,  n} 
and  set  of  edges  E.  Let  d{i,  j)  denote  the  length  of  edge  (i,  j).  We  assume  that 
graph  G{V,  E)  is  complete,  i.e.  any  two  vertexes  are  connected  by  an  edge.  A 
path  on  graph  G{V,  E)  is  defined  as  an  ordered  subset  of  set  V .  The  length  of 
path  P  =  {pi, ...  ,pr}  is  defined  by  the  following  formula: 

r— 1 

L(P)  =  J]d(pi,Pi+i)  (1) 

1=1 

A  path  that  returns  to  the  origin  is  a  cyde.  A  cycle  that  visits  each  vertex  of 
graph  Giy,  E)  once  is  a  tour  on  graph  G{V,  E).  A  tour  (ti  — >  t2  — ^  •  •  •  — ^ 
tr  — »  fi)  is  defined  (and  denoted)  by  an  ordered  set  T  =  {fi,  <2)  •  •  •  >  U}-  The 
length  of  tour  T  is  computed  by  the  following  formula: 

r— 1 

L{T)  =  <iiU,ti+i)  +  d{tr,ti)  (2) 

i=l 

A  cyclic  rotation  of  vertexes  in  a  tour  does  not  change  the  length  of  the 
tour.  A  regular  Traveling  Salesman  Problem  (TSP)  consists  in  finding  the 
shortest  tour  on  graph  G(V,  E): 

Ltsp{N)=  rain  L{T)  (3) 

T-.TcV,\T\=n 

Multi-Traveling  Salesmen  Problem  is  an  extension  of  TSP  to  the  case 
of  several  salesmen.  The  objective  of  the  Multi-Traveling  Salesmen  Problem 
with  m  salesmen  (m-TSP)  is  to  find  a  decomposition  of  graph  G(V,  E)  into  m 
disjoint  complete  subgraphs  minimizing  an  aggregated  objective  resulted  from 
the  particular  solutions  of  TSPs  on  the  subgraphs.  Min-Sum  m-TSP  with  a 
depot  vertex  is  a  prevalent  formulation  of  m-TSP.  According  to  this  formula¬ 
tion,  the  objective  is  to  determine  m  tours  of  the  least  total  length,  so  that  the 
depot  vertex  is  present  in  every  tour,  and  every  other  vertex  of  graph  G(V,  E) 
is  present  in  only  one  tour.  This  variation  of  m-TSP  can  be  transformed  into 
a  regular  TSP  by  introduction  of  artificial  vertexes  [1].  Unfortunately,  the  re¬ 
sultant  TSP  becomes  very  degenerate  in  most  cases.  The  transformation  of 
Min-Sum  m-TSP  with  different  depot  vertexes  appears  to  be  difficult  except 
for  a  special  case  of  only  two  salesmen.  The  transformation  for  m  =  2  is 
provided  in  [5]. 

A  typical  solution  of  Min-Sum  m-TSP  is  highly  irregular.  The  computa¬ 
tional  experience  shows  that  the  distribution  of  assignments  is  not  uniform  in 
the  optimal  solution:  the  difference  in  the  length  of  two  tours  can  be  excessive. 
In  many  contexts  this  can  be  inappropriate:  Giust  [4]  considers  an  example 
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of  gas  distribution  by  a  small  company  that  has  four  delivery  cars.  In  another 
example,  France  et  al.  [2]  describe  a  problem  of  scheduling  n  jobs  to  m  iden¬ 
tical  machines  in  order  to  minimize  the  total  time  when  processing  times  are 
job-sequence  dependent.  The  authors  consider  the  Min-Max  m-TSP  formu¬ 
lation  with  a  depot  vertex.  In  this  formulation  the  objective  is  to  determine 
m  tours  with  the  least  length  of  the  longest  tour;  every  tour  starts  and  ends  at 
the  depot  vertex,  and  every  vertex  of  graph  G{V,  E)  (except  the  depot  vertex) 
is  present  in  only  one  tour.  The  authors  propose  one  heuristic  and  two  exact 
search  schemes  for  the  Min-Max  m-TSP  with  depot  formulation. 

We  consider  Min-Sum  m-TSP  without  a  depot  vertex  (No-Depot  Min- 
Max  m-TSP).  In  this  case,  the  Submarine  Routing  Problem  can  be  considered 
as  a  possible  application  of  this  formulation.  In  this  problem  graph  G(y,  E) 
represents  a  region  that  needs  to  be  monitored.  Each  vertex  of  the  graph  cor¬ 
responds  to  a  specific  geographical  location  in  the  region.  There  is  a  fleet  of 
submarines  available  for  monitoring.  Due  to  the  high  price  of  a  submarine,  the 
fleet  is  very  limited,  and  usually  contains  2-3  submarines.  A  patrol  cycle  of  a 
submarine  is  the  time  needed  to  visit  all  the  assigned  locations  and  to  return 
to  the  origin.  The  objective  of  the  Submarines  Routing  Problem  is  to  assign  a 
specific  route  (tour)  to  each  submarine,  so  that  each  location  is  visited  once  by 
only  one  submarine,  and  the  longest  patrol  cycle  is  minimized.  The  problem 
can  be  formulated  as  follows: 

LmIN-MAX  =  {LTSPiMi)}  (4) 

V  =  (5) 

MiHMj- =  0Vt7^  j  (6) 

In  this  formulation  of  Min-Max  m-TSP  a  depot  vertex  is  not  specified.  The 
submarines  are  assigned  to  their  patrol  routes  for  many  cycles  for  a  period  of 
time  ranging  from  1  to  5  years.  Provision  and  crew  changes  are  provided  in 
several  locations  on  the  patrol  routes  by  special  ships  or  aircrafts.  Usually, 
these  locations  are  chosen  after  the  submarines’  routes  are  specified. 

The  purpose  of  this  study  is  to  study  No-Depot  Min-Max  2-TSP.  We  in¬ 
troduce  a  notion  of  characteristic  function  for  this  class  of  problems.  Using 
constant  graphs  we  study  a  connection  between  No-Depot  Min-Max  2-TSP 
and  a  subclass  of  self-dual  monotonic  Boolean  functions. 

2.  Characteristic  Function  for  No-Depot 
Min-Max  2-TSP 

We  consider  Min-Max  No-Depot  2-TSP  on  complete  undirected  graph 
G{V,  E)  with  non-negative  lengthes  of  edges.  We  consider  a  two-stage  so¬ 
lution  of  the  problem.  At  the  first  stage,  the  partition  of  set  F  =  {1, . . . ,  n} 
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into  two  subsets  Mi  and  M2  is  created.  The  first  set  Mi  is  assigned  to  the  first 
salesman,  and  the  second  set  M2  is  assigned  to  the  second  salesman.  At  the 
second  stage  each  salesman  solves  a  regular  TSP  on  the  subgraph  assigned  to 
him.  We  associate  each  graph  partitioning  {Mi,  M2}  with  a  Boolean  vector 
a  =  {o-i, . . . ,  an}  €  (0, 1}”  in  the  following  way: 


To  if  i  G  Ml 
{1  if  i  G  M2 


(7) 


Further  we  will  identify  each  decomposition  {Mi,  M2}  of  graph  G{V,  E)  by 
Boolean  vector  5  =  {ai , . . . ,  an}  according  to  rule  (7). 

Definition  2.1.  Characteristic  function  f{xi, . . .  ^Xn)  for  No-Depot  Min-Max 
2-TSP  on  graph  G(V,  E)  (or  simply  a  characteristic  function  for  graph  G{V,  E)) 
is  defined  by  the  following  rule: 


f  1,  Lrspi^i)  ^  Ltsp{^2) 
\  0,  Ltsp{M2)  <  Ltsp{M2), 


(8) 


where  a  is  defined  in  (7). 

Boolean  vector  a  =  (ai, . . . ,  an)  is  called  one  if  f{d)j=  1.  One  a  = 
(ai, . . . ,  an)  is  called  a  lower  one  if  there  is  no  other  one  (3  =  (^1, . . . ,  fin) 
so  that  fii  <  ai'ii  =  1, . . . ,  n.  Otherwise  one  a  =  (ai, . . . ,  an)  is  called 
a  generated  one.  Similarly,  Boolean  vector  a  =  (ai, . . . ,  an)  is  called  zero 
if  /(a)  =  0.  Zero  a  =  (ai, . . . ,  an)  is  called  an  upper  zero  if  there  is  no 
other  zero  fi  =  (fii,...,  fin)  so  that  fii  >  ai  for  i  =  1,.. .  ,n.  Otherwise  zero 
a  =  (ai, . . . ,  an)  is  called  a  generated  zero. 

A  graph  is  called  metric  if  all  the  vertexes  in  the  graph  correspond  to  the 
points  in  metric  space,  and  the  lengths  of  the  edges  equal  the  metric  distances 
between  the  corresponding  points.  By  the  definition,  any  metric  graph  is  a 
symmetric  graph,  and  for  any  three  graph  vertexes  ii,  i2,  and  is  the  triangle 
inequality  d{ii,  is)  <  d(ii, *2)  +  d(ii, *2)  is  satisfied. 

Definition  2.2.  Graph  G{V,  E)  is  splittable  if  for  any  three  vertexes  the  trian¬ 
gle  inequality  is  satisfied,  and  for  its  characteristic  function  there  is  no  decom¬ 
position  {Ml,  M2}  of  the  graph,  for  which  Ltsp{M\)  =  LrspiMf). 

(Note,  in  this  definition  it  is  not  necessary  that  the  considered  graph  is  met¬ 
ric) 

Statement  2.1.  The  characteristic  function  of  a  metric  (splittable)  graph  is 
monotonic. 

Indeed,  since  the  triangle  inequality  holds  for  any  three  vertexes  of  a  metric 
(splittable)  graph,  adding  a  new  vertex  to  Mi  does  not  decrease  Ltsp{Mi), 
and  removing  a  vertex  from  M2  does  not  increase  Ltsp{^2)‘ 
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A  Boolean  function  /(5)  is  selMual  if  for  any  Boolean  vector  7+  —  (71 , . . ,  7n) 
and  vector  7"  =  (1  —  71, 1  —  7n) 

/(7+)  +  /(7“)  =  1 

The  next  property  follows  directly  from  definitions  2.1  and  2.2; 

Statement  2.2.  The  characteristic  function  of  a  splittable  graph  is  self-dual. 
Consider  a  non-splittable  metric  graph.  For  this  graph,  the  equality 
Ltsp{Mi)  =  Ltsp{^2) 

holds  for  at  least  one  decomposition  {Mi,  M2}  of  graph  G{V,  E).  The  same 
equation  holds  for  opposite  decomposition  {  M2 ,  Mi } : 

Ltsp{M2)  =  Lrspi^i)- 

Therefore,  the  values  of  the  characteristic  function  are  equal  on  the  vectors 
corresponding  to  decompositions  {Mi,  M2}  and  {M2,  Mi}.  Hence,  the  char¬ 
acteristic  function  of  the  not-splittable  metric  graph  is  not  self-dual.  This  situ¬ 
ation  can  be  easily  avoided  by  small  variation  of  distances  in  the  graph. 

Theorem  2.1.  If  No-Depot  Min-Max  2-TSP  is  considered  for  a  metric  (split- 
table)  graph  G{V,  E),  at  least  one  optimal  solution  of  this  problem  belongs 
to  the  set  of  lower  ones  of  the  corresponding  characteristic  function.  There  is 
another  optimal  solution  that  belongs  to  the  set  of  upper  zeros. 

Proof  Suppose,  that  the  characteristic  function  has  neither  upper  zero,  nor 
lower  one,  that  corresponds  to  an  optimal  solution  of  No-Depot  Min-Max  2- 
TSP  (4)-(6).  Therefore,  every  optimal  solution  of  the  problem  corresponds 
either  to  generaited  zero,  or  to  generated  one.  Because  of  the  symmetry  of  the 
problem,  the  solution  formed  by  the  opposite  decomposition  is  also  optimal.  If 
the  original  solution  corresponds  to  a  zero  (one)  of  the  characteristic  function, 
the  opposite  solution  corresponds  to  a  one  (zero)  of  this  function.  Consider 
optimal  solution  a®  =  (a?, . . . ,  a°)  that  corresponds  to  a  zero  of  characteristic 
function  /(a): 

0  _  /  0  if  i  €  Ml 

\  1  ifi€M2  ’ 

Ml  U  M2  =  V,  Ml  n  M2  =  0. 

According  to  the  assumption,  a®  is  a  generated  zero.  Therefore,  there  exists  an 
upper  zero  cT*  that  exceeds  a®  in  several  components.  So,  there  is  a  non  empty 
set  of  vertexes  H,H  C  M2,  so  that: 


(  0  ifieMiUi? 

\  1  ifi€M2\H  ’ 
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According  to  Definition  2. 1  for  zeros  a*  and  a®  the  following  inequalities  are 
satisfied: 

Ltsp{M\)  ^  Ltsp{M2),  (9) 

Lj'Sp{M\  U  iy)  ^  Lxsp{^2  \  H').  (10) 

Since,  the  graph  is  metric  (splittable),  according  to  Statement  2.1  the  following 
conditions  are  valid: 

L7’5p(Mi)  <  Ltsp(-^i  U -ff))  (1^) 

Ltsp{^2  \  H)  <  Ltsp{^2)-  (12) 

Hence,  we  immediately  obtain: 

min{Lr5p(Mi  U  H),  Ltsp{M2  \  H)}  <  min{Lpsp(M),  i'TSp(-W2)}. 

(13) 

Since,  by  the  proposition,  the  value  of  expression 

min{LT5p(Mi),  Ltsp{M2)} 

is  minimal  for  all  bi-partitions  of  set  V,  the  following  equality  is  satisfied; 

min{L7’5p(Mi  U  H),  Ltsp{M2  \  H)} 

=  min{Lrsp(-W^i),  Lt5p(M2)}- 

Therefore,  for  characteristic  function  f{a)  there  exists  upper  zero  a*  corre¬ 
sponding  to  an  optimal  solution  of  No-Depot  Min-Max  2-TSP  problem  (4)-(6). 
Opposite  vector  a**  =  1  —  a*  is  a  lower  one  of  the  characteristic  function;  this 
vector  corresponds  to  another  optimal  solution  of  the  problem. 

□ 

According  to  Theorem  2.1  each  No-Depot  Min-Max  2-TSP  problem  has  a 
corresponding  self-dual  monotonic  Boolean  function.  In  the  following  sections 
we  consider  the  reverse  question: 

What  Boolean  function  has  a  corresponding  No-Depot  Min-Max  2-TSP,  for 
which  this  function  is  characteristic! 

Below  we  demonstrate  that  for  every  threshold  self-dual  monotonic  func¬ 
tion  it  is  possible  to  find  a  No-Depot  Min-Max  2-TSP,  for  which  the  consid¬ 
ered  function  is  characteristic.  For  the  case  discussed  in  the  next  sections  the 
developed  graph  is  splittable. 

3.  Threshold  Characteristic  Function 

In  this  section  we  consider  Boolean  functions  defined  on  set  {—1, 1}”.  This 
definition  of  Boolean  function  is  different  from  the  standard  one  when  the  func¬ 
tion  is  defined  on  {0, 1}".  We  use  this  format  to  emphasize  specific  properties 
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of  the  considered  problem.  Moreover,  all  the  necessary  properties  of  standard 
Boolean  function  remain  valid  for  this  class  of  functions.  We  identify  every 
Boolean  vector  0  =  {0i, ,  Pn)  by  vector  \\  Pi,.  ..,Pn  II  In  linear  space  L„. 

A  set  of  vectors  on  which  /(/?)  =  0  is  called  a  set  of  zeros  of  function 
f{P),  and  is  denoted  by  f~H0).  Set  of  ones  /”Hl)  ‘s  defined  similarly. 
Boolean  function  is  called  threshold  [6],  if  there  exists  a  set  of  real  numbers 
xi,X2, ...  ,Xn,c,  so  that  linear  inequality 

XiPi  +  X2P2  +  •  •  •  +  XnPn  <  C 

holds  for  Boolean  vector  P  =  {Pi,.  ■  ■  ,Pn)  tfffiP)  =  0-  Thus,  for  threshold 
Boolean  function  f{0)  sets  /-^O)  and  are  separated  by  hyperplane 

XiPi  +  X2P2  +  •  •  •  +  XnPn  —  C) 

which  is  called  a  threshold  hyperplane  for  threshold  Boolean  function  /(/?). 
A  hyperplane  is  called  a  central  hyperplane  if  the  coordinate  origin  belongs  to 
this  hyperplane. 

The  following  statement  makes  a  connection  between  self-dual  threshold 
Boolean  functions  and  central  hyperplanes. 

Statement  3.1.  Threshold  Boolean  Junction  f  {0)  is  self-dual  if  and  only  if 
there  exists  a  central  threshold  hyperplane  for  this  function. 

Indeed,  presence  of  a  central  threshold  hyperplane  implies  self-duality  of 
f{P).  Suppose  now  that  threshold  hyperplane  a;i/3i+X2/32+-  ..-^-XnPn  =  cof 
self-dual  Boolean  function  f{0)  is  not  passing  through  the  coordinate  origin. 
Because  of  the  self-duality,  sets  f~^  {0)  and  are  symmetric  to  each 

other.  Consequently,  hyperplane  si^i  +  X2P2  -!-•••  +  XnPn  =  ~c,  which  is 
symmetric  to  the  original  threshold  hyperplane,  is  also  a  threshold  hyperplane. 
Since  the  considered  hyperplanes  are  threshold  at  the  same  time,  the  following 
conditions  are  valid: 

if  f  {P)  =  1  then  xiPi  X2P2  -H  •  ■  •  -b  XnPn  ^  c; 

if  f{0)  =  0  then  xiPi  -\-  X2P2  -!-•••  +  XnPn  < 

Consequently,  central  hyperplane 

XiPi  -f  X2P2  +  •  •  •  +  XnPn  =  0 

is  a  threshold  hyperplane  for  Boolean  function  f{0).  Therefore,  we  obtained  a 
central  threshold  hyperplane  that  divides  linear  space  Ln  into  two  areas  Aq  = 
{0  :  XiPi+X2P2+.  •  .+XnPn  <  0}  and  Ai  =  {0  :  XiPi+X2P2  +  .  •  .+XnPn  > 
0}. 
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Hence,  the  next  statement  follows  directly: 

Statement  3.2.  /~^(0)  C  Aq  and  C  Ai. 

Indeed,  because  of  self-duality,  neither  zero,  nor  one  of  Boolean  function 
f0)  belongs  to  the  central  threshold  hyperplane.  Therefore,  all  zeros  of 
Boolean  function  f0)  belong  to  area  Aq,  whereas  all  ones  of  the  function 
belong  to  area  A\. 

Definition  3.1.  Zero  (/3i, . . . ,  A-i,  -1,  A+i)  ■■•,Pn)  of  monotonic  Boolean 
function  f0)  is  the  frontier  zero,  if  (/0i, . . . ,  A-ii  +!>  A-fij  •  •  •  > Pn)  is 
a  one  of  function  f0).  One  (/?i, . . . ,  A-i,  H-l,  Pi+i,  ■■■,Pn)  of  monotonic 
Booleanfunction  f{P)  is  the  frontier  one,  ,  j3i-i,—l,  Pi+i, . . . ,  Pn) 

is  a  zero  of  function  f{P). 

The  next  statement  follows  immediately: 

Statement  3.3.  For  a  monotonic  self-dual  Boolean  function  there  exists  at  least 
one  frontier  zero,  and  at  least  one  frontier  one. 

Definition  3.2.  Function  f0)  is  essentially  independent  of  the  variable,  if 
for  any  Pi,...,  Pi-i,  Pi+i,  ...,Pn^  {-1,  +1} 

f{.Pli  •  •  •  >  Pi—ti  1)  Pi+lt  •  •  •  )  Pn)  ~  f{.Plt  •  •  •  1  Pi—ti  "I”!}  /^i+1)  •  •  •  >  Pn)j 

(15) 

otherwise,  function  /(/?)  is  essentially  dependent  on  the  variable, 

Statement  3.4.  Monotonic  Booleanfunction  f{P)  is  essentially  dependent  on 
the  variable  if  and  only  if  this  function  has  the  frontier  zero  and  the  z**' 
frontier  one. 

This  statement  immediately  follows  from  Definitions  3.1  and  3.2. 

Theorem  3.1.  For  any  threshold  self-dual  Booleanfunction  f{P),  essentially 
dependent  on  all  its  variables  and  defined  by  threshold  hyperplane  xifii  + 
X2P2  +  • . .  +  XnPn  =  c  nontrivial  vector  ||  oi, . . . ,  o„  ||,  ai, . . . ,  On  >0 
belongs  to  area  Ai  =  {P  '.  xiPi  -f  X2P2  -!-•••  +  XnPn  >  0},  whereas  the 
opposite  vector  ||  — oi, . . . ,  —an  ||  belongs  to  area  Aq  =  {P  :  xiPi  ■+■  X2P2  + 

. . .  -t-  XnPn  <  0}. 

Proof.  At  first,  we  show  the  validity  of  the  theorem  for  vector  ||  1, 0, . . . ,  0  ||. 
Denote  R=  {P  :  xiPi  +  X2P2  +  • . .  +  XnPn  =  0}.  Accordingjo  Statement 
3.1,  is  a  central  threshold  hyperplane  for  Boolean  function  /(/?). 

Suppose  that  vector  ||  1, 0, . . . ,  0  ||  does  not  belong  to  Ai.  Therefore, 

II  1,0, . . . ,0  ||€  Ao U i?. 
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If  II  1, 0, . . . ,  0  II €  R  then  xi  =  0,  and  hyperplane  X202  +  •  XnPn  =  0 
is  a  threshold  hyperplane  of  function  f0).  Hence,  function  / (0)  does  not  de¬ 
pend  on  variable  (3%,  that  contradicts  the  condition  of  the  essential  dependency. 
Therefore 

II  l,0,...,0||€^o. 

Since  function  /(/?)  is  essentially  dependent  on  all  its  variables,  due  to 
Statement  3.4,  this  function  has  the  1-st  frontier  zero  ||  —1,  ^2)  •  •  •  >  /^n  11^  -^o- 
Since  Aq  is  a  conical  set,  and  ||  1, 0, . . . ,  0  ||  and  ||  — 1,  y02,  •  •  •  >  Pn  ||  belong 
to  that  set, 

(II  1,0,...,0||  +2  II  -l,^2,...,^n||)e^0. 

By  Definition  3.1, 

II  11=11  1,0,... ,0  II  -h2  II  -l,/32,...,/3n  || 

is  a  frontier  one,  and,  consequently,  belongs  to  set  Ai.  Therefore, 

Ai  U  j4o  ^  0. 

Thus,  we  obtain  a  contradiction  of  what  we  have  assumed.  Therefore, 

II  l,0,...,0||€^i. 

Similarly,  the  statement  can  be  proved  for  any  other  unit  vector 
II  0,...,0,1,0,...,0||  . 

Since  Ai  is  a  conical  set,  any  nontrivial  linear  conical  combination  (combina¬ 
tion  with  non-negative  coefficients)  of  its  elements  belongs  to  this  set.  There¬ 
fore,  for  any  nontrivial  combination  of  nonnegative  numbers  ai, . . . ,  Cn, 

II  — Ol,  •  •  ■ ,  ||€  Aq 


and 


II  ai, . . .  (On  ||€  A\. 


□ 


The  following  statement  is  an  immediate  inference  of  the  theorem: 

Statement  3.5.  If  threshold  self-dual  Boolean  function  f  0)  is  essentially  de¬ 
pendent  on  all  its  variables,  and  xifii  -|-  X2P2  +  •  ■  ■ + XnPn  =  c  is  the  threshold 
hyperplane  of  the  function  f0),  then  Xi  >  for  any  i  =  1, . . . ,  n. 

Indeed,  according  to  Theorem  3.1 ,  the  vectors  of  form  ||  0, . . . ,  0, 1, 0, . . . ,  0  || 
satisfy  inequality  xiPi  -I-  X2/h  +  .  •  •  +  XnPn  >  0,  that  proves  the  statement. 
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4.  Constant  Graphs 

Definition  4.1.  Symmetric  graph  G{V,  E)  is  a  constant  graph  if  there  exists  a 
set  of  numbers  B  =  {6i  so  that  for  every  pair  of  vertexes  i  and  j  of  graph 

G{V,  E),  the  length  of  a  connecting  edge  is  d{i,  j)  =  6i  +  bj.  Set  B  = 
is  a  defining  set  for  constant  graph  G{V,  E). 

The  TSP  considered  on  a  constant  graph  belongs  to  the  class  of  Constant 
Discrete  Programming  Problems,  which  are  studied,  for  example,  in  [3]. 

Statement  4.1.  For  constant  symmetric  graph  G(V,  E)  the  length  of  a  cycle 
does  not  depend  on  the  order  of  the  visits  to  the  vertexes  of  the  cycle. 

Indeed,  the  length  of  cycle  c  =  {ji  ^  j2  ^  3k  ji)  is  determined 
by 

k-l 

He)  =  YlH3h3i+i)  +  dUkJi)  (16) 

i=l 

k-l  k-l 

—  +  ^ji+i )  +  +  bji  =  2  ^  =  2  ^  bj^.  (17) 

z=l  i=l 

Moreover,  it  can  be  proved  that  if  for  a  complete  graph  the  length  of  a  cycle 
does  not  depend  on  the  order  of  the  visits  to  the  vertexes  of  the  cycle,  the  graph 
is  constant  [3].  As  a  direct  consequence  of  Statement  4.1,  TSP  is  a  trivial 
problem  for  a  constant  graph.  Indeed,  every  tour  in  that  graph  has  the  same 
length. 

Statement  4.2.  For  constant  graph  G{V,  E)  with  non-negative  elements  of 
defining  set  B  =  {6t}iZy,  the  triangle  inequality  holds  for  all  triplets  of  ver¬ 
texes  of  that  graph. 

Indeed, 

d{i,j)  +  d{j,  k)  =  (bi  +  bj)  +  (bj  +  bk) 

—  (pi  "b  bk)  ~\~  ^bj  —  d{i,  k)  +  2bj  ^  d(i,  k) 


5.  Interpretation  of  Threshold  Self-Dual 
Monotonic  Boolean  Functions 

Consider  symmetric  constant  graph  G{V,  E)  with  nonnegative  elements  in 
its  defining  set  B  =  •  As  it  has  been  shown  in  the  previous  section,  the 

triangle  inequality  holds  for  any  triplet  of  the  graph  vertexes,  and  the  length  of 
a  cycle  does  not  depend  on  the  order  of  the  visits  to  the  vertexes. 

Assuming  that  graph  G{V,  E)  is  splittable,  consider  No-Depot  Min-Max 
2-TSP  assigned  to  this  graph.  The  characteristic  function  of  constant  graph 
G(V,  E)  is  defined  by  the  following  mle: 
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where 


if  2  2  YjieM2 

otherwise  ’ 


„  _  /  0  if  i  €  Ml 
\  1  ifieMa  ’ 

Ml  U  M2  =  V,  Ml  n  M2  =  0. 


(18) 

(19) 


Using  the  vector  notations  introduced  above  and  an  assumption  that  graph 
G{V,  E)  is  splittable,  the  definition  of  function  f{0)  can  be  rewritten  in  the 
following  way: 


f0)  = 


1  if;^x>0 
0  if^x<Q 


(20) 


where 

P  —  11^1)  •  •  • )  Ai||>  ®  ~  II®1>  ■  •  • )  ^n\\' 

Since  characteristic  function  /(/3)  is  monotonic,  the  determination  of  all  its 
lower  ones  is  sufficient  to  completely  describe  the  function.  Let  Fi  denote  the 
set  of  lower  ones ,  and  Fq  denote  the  set  of  upper  zeros  of  function  /  (/?) .  Due  to 
the  self-duality  of  function  f{/3)  the  following  statement  is  valid  for  splittable 
graph  G{V,E): 

Statement  5.1.  ^  €  Fi  -0  e  Fq,  that  is,  the  vector  opposite  to  a  lower 
one  is  an  upper  zero,  and  vice  versa. 

Definition  5.1.  Afundamentcd  matrix  of  ones  of  monotonic  Boolean  function 
f0)  is  matrix  Di,  constructed  by  the  following  rules:  the  first  n  rows  form 
an  identity  matrix  n  x  n;  the  next  rows  are  formed  by  all  the  vectors  from 
Fi  arranged  in  binary  increasing  order.  A  fundamental  matrix  of  zeros  of 
monotonic  Boolean  function  f0)  is  matrix  Dq,  constructed  by  the  following 
rules:  the  first  n  rows  form  a  negative  identity  matrix  n  x  n;  the  next  rows  are 
formed  by  all  the  vectors  from  Fq  arranged  in  binary  decreasing  order. 

Due  to  self-duality  of  function  f0),  D\  —  —Dq. 

Statement  5.2.  Any  one  0i  of  monotonic  Boolean  function  f0)  can  be  rep¬ 
resented  as  pi  =  116i, . . . ,  6m||  X  Di;  any  zero  fio  of  monotonic  Boolean  func¬ 
tion  f0)  can  be  represented  as  fio  =  ||(>ij  •  •  •  i  ^m||  x  Dq,  where  bi  €  {0, 2}, 
i  =  1, . . . ,  n,  and  6i  €  {0, 1},  i  =  n  -M, . . .  ,m,  and  m  is  the  number  of  rows 
in  fundamental  matrix  D\  (Dq). 

Indeed,  one  of  function  f0)  is  generated  by  some  lower  one  Sup¬ 
pose  that  Pi  is  represented  by  the  fe***  row  of  fundamental  matrix  Di.  Then 

A  =  ||6i,...,&n,0,...,0,6fc  =  l,0,...,0||  xDi,  (21) 
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where 

^  _  f  2  if  01  differs  from  0^  in  the  i**’  component  ^22) 

*  1  0  if  /5i  conicides  with  0i  in  the  component 

Similarly,  it  can  be  shown  for  zero  0o  of  function  f{0). 

For  the  rest  of  the  section,  notation  Dxx>0(Dxx<0)  means  that 
resultant  vector  v  =  D  x  xhas  only  positive  (negative)  components. 

Theorem  5.1.  Self-dual  monotonic  Boolean  function  /(/?)  essentially  depen¬ 
dent  on  all  its  variables  is  an  threshold  if  and  only  if  there  exists  a  vector  x, 
for  which  Di  x  x  >  0,  where  Di  is  a  fundamental  matrix  of  ones  of  function 

f0)- 

Proof  Existence:  Due  to  Statement  3.1  self-dual  threshold  function  f{0)  has 
central  threshold  hyperplane  {0,  x)  =  0,  and  for  any  one  /?i  of  function  /(/3) 

{01,  x)  >  0. 

According  to  Statement  3.5,  all  components  of  vector  x  are  positive.  Since  the 
rows  of  matrix  D\  are  either  rows  of  the  unit  matrix,  or  lower  ones  of  function 
f{0),  the  following  inequality  is  satisfied: 

D\xx>f). 

Self-Duality:  Suppose  that  there  exists  vector  x,  so  th^  Dix  x>  0.  Con¬ 
sider  any  one  0\  of  function  f{0).  Due  to  Statement  5.2,  Pi  =  \\bi, . . .  ,bm\\^ 
Di,  where  bi  €  {0, 2},  i  =  1, . . . ,  n,  and  bi  €  {0, 1},  i  =  n  -t- 1, .  -  • ,  m,. 
Denote  v  -  Di  x  3.  Since  Ui  >  0  and  6*  >  0  for  any  i  =  1, . . . ,  m,  and  there 
exist  i*,  for  which  bi*  >  0,  then 

(/3i,x)  =  (6  X  Di,x)  =  (6,  (£>i  X  xf)  > 

>  bi*  X  min  >  0. 
i 

Hence,  for  one  0i  of  function  f{0) 

{01,3)  >Q. 

Since  function  f{0)  is  self-dual,  for  any  zero  0  of  function  /  (/5) 

(y^,x)<0. 

Therefore,  function  f{0)  is  a  threshold  one  with  hyperplane  {0, 3)  =0.  □ 
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Theorem  5.2.  For  any  threshold  self-dual  monotonic  Boolean  function  /(^) 
essentially  dependent  on  all  its  variables  there  exists  a  graph,  for  which  func¬ 
tion  f{P)  is  characteristic. 

Proof.  Consider  threshold  self-dual  monotonic  Boolean  function  / (/?),  which 
is  essentially  dependent  on  all  its  variables.  Let  D\  be  a  fundamental  matrix 
of  ones  of  function  f0).  Due  to  Theorem  5.1,  there  exists  vector  x  for  which 
Dixx  >  0.  Since  the  first  n  rows  of  matrix  Di  form  an  identity  matrix,  the 
first  n  inequalities  in  Di  x  x  >  0  set  x  >  0.  The  rest  of  Ae  rows  of  matrix 
D\,  by  definition,  correspond  to  lower  ones  of  function  /(/0).  Define  constant 
graph  G{y,  E)  by  defining  set  B  =  {bi  =  Xi,  i  =  1, ... ,  n}.  The  length  of 
edge  (ij)  is  d{i,j)  =  bj  =  Xi-\-  xj.  Consider  No-Depot  Min-Max  2- 

TSP  on  this  graph.  Let  f0)  be  a  characteristic  function  of  the  graph.  Due  to 
Statements  2.1  and  2.2  function  f0)  is  monotonic  and  self-dual.  Therefore, 
the  following  condition  satisfies  each  one  Pi  (and  each  zero  Po  =  —Pi)  of 
function  f{P): 

0ux)  >  0. 

Since  Z>i  x  s  >  0  and  all  the  lower  ones  of  function  f{P)  are  among  the  rows 
of  matrix  Di,  every  lower  one  of  Boolean  function  f{P)  is  a  one  of  function 

f{P).  Because  of  monotonicity  and  self-duality  of  functions  f{p)  and  f{P), 
and  the  fact,  that  any  monotonic  Boolean  function  is  completely  defined  by  the 
set  of  lower  ones,  the  equation 

m = R0) 

is  valid  for  any  binary  0.  □ 

Note,  that  in  the  proof,  constructed  graph  G{V,  E)  is  splittable. 

Consider  the  case  when  Boolean  threshold  function  f{P)  is  independent  of 
some  of  its  variables.  The  following  statement  is  valid  for  this  function: 

Statement  5.3.  If  threshold  function  f{P)  is  independent  of  the  i^^^variable, 
and  {P,x)  =  c  is  a  threshold  hyperplane  for  this  function,  then  {P,y)  =  c, 
where 

_  j  0  j  =  i 
\  Xj  j  i 

is  also  a  threshold  hyperplane  of  function  f{P) 

Indeed,  consider  a  one  pi  of  function  f{P).  It  satisfies  (Pux)  >  c.  Since 
function  f{P)  is  independent  of  the  variable,  there  exists  a  one  Pi  of  the 


Properties  of  No-Depot  Min-Max  2-Traveling-Salesmen  Problem 


97 


function  that  differs  from  p  by  only  the  variable.  Therefore,  {0i  ,x)  >  c. 
Adding  inequalities  (A>®)  ^  c  and  ,x)  >  c,  jve  get  that  ^  satisfies 
{PiiV)  >  c.  Similarly,  any  zero  0q  of  function  f0)  satisfies  {0o,y)  <  c. 
Therefore,  {0,  ^  =  c  is  a  threshold  hyperplane  of  function  /(/3).  The  reverse 
statement  appears  straightaway: 

Statement  5.4.  If  threshold  Boolean  function  f{0)  has  threshold  hyperplane 
0,  x)  =  c,  and  Xi  =  0,  then  function  f0)  is  independent  of  the  variable. 

Definition  5.2.  Central  threshold  hyperplane  0,  x)  =  Oof  threshold  Boolean 
function  f0)  is  a  reduced  hyperplane,  if  for  every  variable,  of  which  function 
f0)  is  independent,  the  corresponding  component  of  vector  x  is  zero. 

A  central  threshold  hyperplane  of  a  Boolean  function  essentially  dependent 
on  all  its  variables  is  a  simple  case  of  a  reduced  hyperplane.  Due  to  Statements 
3.1  and  5.3  and  the  definition  of  central  threshold  hyperplane,  the  following 
statement  is  valid: 

Statement  5.5.  For  any  threshold  self-dual  Boolean  function  there  exists  a 
reduced  central  threshold  hyperplane. 

We  are  ready  to  prove  the  main  result  of  the  chapter. 

Theorem  5.3.  For  any  threshold  self -dual  monotonic  Boolean  function  f{P) 
there  exists  a  graph,  for  which  function  f{(3)  is  characteristic. 

Proof.  The  case  with  function  f0)  essentially  dependent  on  all  its  variables 
has  been  considered  in  Theorem  5.2.  Now  we  consider  the  case  when  function 
f0)  is  independent  of  some  of  its  variables.  Let  Kindep  =  {1,  • . . ,  fc}  be  a 
set  of  indexes  of  variables,  on  which  the  function  /(/3)  does  not  depend,  and 
Kdep  =  {A:  +  1, . .  • ,  n}  be  a  set  of  indexes  of  the  variables,  on  which  function 
f0)  essentially  depends.  According  to  Statement  5.5,  there  exists  reduced 
central  threshold  hyperplane  0,^  =  0  for  function  f{P).  Define  constant 
graph  G{V,  E)  by  defining  set  B  =  {bj  =  Vi  i  =  1,  •  •  •  >  ’T'}-  The  length  of 
edge(i,j)is 

d{i,  j)  =  bi  +  bj  =  Vi  +  Vj. 

For  any  variable,  on  which  function  f0)  does  not  depend,  the  corresponding 
component  of  y  is  zero.  In  compliance  with  the  proof  of  Theorem  5.2,  Boolean 
function  f0)  is  a  characteristic  function  for  No-Depot  Min-Max  2-TSP  con¬ 
sidered  on  graph  G(V,£?).  D 

The  same  as  in  the  proof  of  Theorem  5.2,  the  obtained  graph  is  splittable. 
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6.  Conclusion 

The  Multiple  Traveling  Salesmen  Problem  has  many  variations.  In  this  work 
we  studied  a  No-Depot  Min-Max  formulation  of  2-Traveling  Salesmen  Prob¬ 
lem.  We  introduced  a  characteristic  function  for  this  class  of  problems.  This 
Boolean  function  is  monotonic  and  self-dual  for  complete  graphs  with  metric 
distances.  For  an  arbitrary  monotonic  threshold  self-dual  Boolean  function  we 
have  proven  existence  of  a  No-Depot  Min-Max  2-Traveling  Salesmen  Prob¬ 
lem,  for  which  this  function  is  characteristic. 
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